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Abstract. Motivated by recent experiments with confined binary liquid mixtures 
near demixing, we study the universal critical properties of a system, which belongs 
^ ■ to the Ising universality class, in the film geometry. We employ periodic boundary 

conditions in the two lateral directions and fixed boundary conditions on the two 
confining surfaces, such that one of them has a spatially homogeneous adsorption 
preference while the other one exhibits a laterally alternating adsorption preference, 
resembling locally a single chemical step. By means of Monte Carlo simulations 
of an improved Hamiltonian, so that the leading scaling corrections are suppressed, 
^J , numerical integration, and finite-size scaling analysis we determine the critical Casimir 

force and its universal scaling function for various values of the aspect ratio of the 
£Nj ■ film. In the limit of a vanishing aspect ratio the critical Casimir force of this system 

^ | reduces to the mean value of the critical Casimir force for laterally homogeneous -\ — \- 

and H — boundary conditions, corresponding to the surface spins on the two surfaces 
Q^ , being fixed to equal and opposite values, respectively. We show that the universal 

C*") ' scaling function of the critical Casimir force for small but finite aspect ratios displays 

[~^. ■ a linear dependence on the aspect ratio which is solely due to the presence of the 

^D \ lateral inhomogeneity. We also analyze the order-parameter profiles at criticality and 

their universal scaling function which allows us to probe theoretical predictions and to 
compare with experimental data. 






m 



o 



- 1—^ 

X 

^ . PACS numbers: 05.70.Jk, 64.60. an, 68.15. +e, 05.50.+q, 05.10.Ln 

AMS classification scheme numbers: 82B27, 82B20, 76A20, 82B80 



Critical Casimir forces and adsorption profiles 2 

1. Introduction 

If a fluctuating field is confined between surfaces, effective forces arise between them. 
In a system close to a second-order phase transition, the order parameter develops such 
long-ranged fluctuations. The resulting effective force is known as the critical Casimir 
force. This phenomena, first predicted by Fisher and de Gennes [1] is the analogue of the 
Casimir effect in quantum electrodynamics [2]. Reference [3] provides a recent review 
which illustrates analogies as well as differences between the two effects and guides the 
reader towards further reviews of the subject and the pertinent original literature. 

The critical Casimir force is characterized by a universal scaling function, which 
is determined by the bulk and surface universality classes (UC) [HE] of the confined 
system. It is independent of microscopic details of the system, and it depends only on 
a few global and general properties, such as the spatial dimension d, the number of 
components of the order parameter, the shape of the confinement, and the boundary 
conditions (b.c.) there [EHE]. We are solely interested in isotropic systems because the 
critical Casimir force is only active in fluid systems which allow the ordering degrees of 
freedom to move in and out of the system; this is not the case for magnets for which the 
critical Casimir force is not a measureable quantity. (For a discussion of the potential 
relevance of lattice anisotropies with respect to universality JUITO] see Ref. [TT].) 

In recent years the critical Casimir effect has attracted numerous experimental 
[T2TfT8] and even more theoretical investigations. Critical Casimir forces can be obtained 
indirectly by studying wetting layers of fluids close to a critical end point [191120] . In 
this context, 4 He wetting films close to the onset of superfluidity [12] and wetting 
layers of classical [13] and quantum [H] binary liquid mixtures have been studied 
experimentally. Only recently direct measurements of the critical Casimir force have 
been reported [T5HT8] by monitoring individual colloidal particles immersed into a binary 
liquid mixture close to its critical demixing point and exposed to a planar wall. 

Until recently theoretical investigations of the Casimir force used, to a large 
extent, field-theoretical methods. The geometry studied most is the film geometry 
in the presence of symmetry-conserving boundary conditions, which has been studied 
within mean field theory, the e— expansion, and the perturbative expansion at fixed 
dimension d = 3 [T01I2TH35] . real-space renormalization-group theory [19], and within 
the large- N expansion [26 | l2T tl36] - Hd] . Film geometries with symmetry-breaking surface 
fields have been considered within mean field theory [23lll2j|l3] and by using a local 
density functional method [44] . In dimension two exact results are available for periodic 
b.c. [45J, symmetry-conserving b.c. [46] and for symmetry-breaking b.c. [47]. The film 
geometry in the presence of inhomogeneous b.c. (i.e., with patterned substrates) has 
been investigated within mean field theory [48] and within Gaussian approximation 
[4"9] . Curved boundaries, which are of particular experimental interest, pose additional 
difficulties. Spherical geometry has been studied within the e— expansion at criticality 
[50] and, off criticality, within mean field theory, the small-sphere expansion and the 
Derjaguin approximation [5T1453] . The critical Casimir force has been studied also for 
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geometrically structured confinements [SUES] and for ellipsoidal colloids [55JE7] within 
mean field theory. 

Early numerical simulations for the critical Casimir force have been employed in 
Ref. [23] for laterally homogeneous b.c. Quantitatively reliable determinations of critical 
Casimir forces for laterally homogeneous b.c. have been obtained recently by means of 
Monte Carlo simulations. In this context, the XY UC [2? t I58H63] . which describes 
the critical properties of the superffuid phase transition in 4 He, as well as the Ising 
UC J27JEH1E0JEI], which describes the demixing transition in a binary liquid mixture, 
have been investigated and the critical Casimir force has been determined. 

Experiments with binary liquid mixtures have also been used to study critical 
Casimir forces acting on a colloid close to a chemically structured substrate [HI HE], 
creating a laterally varying adsorption preference. Such a system has been investigated 
theoretically for the film geometry within mean field theory [IE] and for a sphere near a 
planar wall within the Derjaguin approximation |53j. Motivated by these experimental 
results, and in the view of the lack of bona fide theoretical data in d = 3 for laterally 
inhomogeneous critical systems, here we present a Monte Carlo study of a three- 
dimensional lattice model in the film geometry, representing the Ising UC. We employ 
periodic boundary conditions in the two lateral directions and fixed boundary conditions 
on the two confining surfaces. The Ising spins on the upper surface are fixed to +1. 
The lower surface is divided into two halves, one with spins fixed to —1 and the other 
with spins fixed to +1, such that the system remains translationally invariant in one 
lateral direction. The lattice constant of the simple cubic lattice is set to 1. Locally this 
mimics a single chemical step as the simplest building element for lateral heterogeneities. 
However, the presence of lateral periodic b.c. in the direction normal to the chemical 
step generates an additional chemical step at the lateral boundaries, resulting in a pair of 
individual chemical steps. In a system of finite lateral size, the presence of these chemical 
steps gives rise to a line contribution to the free energy. We note that in a system with 
a finite lateral size these line contributions cannot be due to a single chemical step 
only. For instance, in the presence of open lateral b.c, besides the single chemical step 
there is an additional line contribution to the free energy due to the free edges. Our 
choice of laterally periodic b.c. corresponds to the simplest implementation of a model 
with a laterally varying adsorption preference. Furthermore, here we consider also the 
extrapolation to an infinite lateral extent; in this limit, the lateral b.c. are irrelevant 
and we recover locally a single chemical step. For this system we determine the critical 
Casimir force and the order-parameter profiles. Certain preliminary results for this kind 
of system have been reported in Ref. [65J. 

Taking corrections-to-scaling into account is important in order to be able to 
extrapolate data for systems of finite size L to the thermodynamic limit L — > oo. 
In particular, in the standard three-dimensional Ising model, scaling corrections are 
proportional to L~ u , with uj = 0.832(6) |66j. The presence of non-periodic boundary 
conditions such as in the film geometry gives rise to additional scaling corrections, the 
leading ones being proportional to L^ 1 , which are numerically difficult to disentangle 
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from the previous ones. In order to avoid the simultaneous presence of these competing 
corrections, we have studied a so-called improved model [67], for which the leading 
scaling corrections oc L~^ are suppressed for all observables so that the correction oc L^ 1 
becomes the leading one. 

This paper is organized as follows. In section |2] we recall the finite-size scaling 
theory for a system with non-periodic boundary conditions, paying special attention to 
the corrections-to-scaling. In doing so we consider the film geometry described above, as 
well as more general b.c. In section[3]we describe the method we use, which is based on 
a combination of Monte Carlo simulations and numerical integration. Our results for the 
Casimir force amplitude at the critical temperature and for the corresponding universal 
scaling function are presented in section H] and section |5l respectively. In section [6] we 
analyze the order-parameter profiles. We summarize our results in section [TJ In order 
to estimate non-universal amplitudes, the analysis of the high-temperature series of the 
improved model is presented in Appendix A In Appendix B we give some technical 



details concerning our Monte Carlo simulations. 

2. Finite-size scaling 

2.1. General properties 

We consider a three-dimensional system in a box of size L x L\\ x L\\, which in 
the thermodynamic limit exhibits a second-order phase transition at the temperature 
T = T c . In this subsection we provide those finite-size scaling (FSS) properties which 
hold independently of the choice of the boundary conditions (b.c), leaving the discussion 
of their influence to section [2T2l Away from criticality (i.e., L, L\\ S> £ ~ \t\~ v , with £ as 
the bulk correlation length, t = (T—T c )/T c , and v the thermal critical exponent), the free 
energy density T per k B T of the system (i.e., the free energy divided by Ll2k B T) can 
be expanded into specific geometric contributions, corresponding to the bulk, surfaces, 
lines, and corners [6H]- All these terms but the bulk one depend on the b.c. (see the 
discussion in the following subsection). 

In the critical region and in the presence of an external bulk field H, the free energy 
density can be decomposed into a singular contribution and a non-singular background 
term: 

T{% H, L, L,|) = J^ (s) (t, H, L, L,|) + ^ ns \t, H, L, L {] ). (1) 



As we shall discuss in section |2T2| analogously to the expansion away from criticality 
the non-singular background ^^(t, H, L, L\\) can be decomposed into specific geometric 
terms. According to renormalization-group (RG) theory [69], neglecting corrections to 
scaling, in spatial dimension d the singular part of the free energy density obeys the 
following scaling property: 

^ s \t,H,L,L {l ) = ±- d f{utV>\u h LV\p), (2) 
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where y t and y^ are the RG-dimensions of the non-linear scaling fields u t and Uh 
associated with the deviations from the critical temperature (u t ~ t) and with the 
external field (uh ~ H), respectively, and p = L/L\\ is the aspect-ratio. The RG- 
dimensions y t and yh are related to standard critical exponents by y t = 1/u and 
yh = (d + 2 — rj)/2. The scaling fields are analytical functions of the parameters of 
the Hamiltonian [691170] which can be expanded as 

u t = a t + o(t,H), (3) 

and 

u h = b H + o(t,H), (4) 

where a^ and &o are non-universal constants. Close to the bulk critical point, the bulk 
correlation length £ varies as 

£(t->0 ± )=£ ± |t|^, (5) 

where the signs ± apply for T ^ T c and £^ are non-universal amplitudes forming the 
universal ratio £o"/£o~ = 1-896(10) [71]. Here and in the following £ is the so-called 
exponential (or true) correlation length, which governs the exponential decay of the 



two-point correlation function (see Appendix A ). Furthermore, one can fix ao as 

ao=(tir 1/U , (6) 

so that the scaling variable is given by 

— LK - t {0=^<i)\i) i - (7) 

The scaling function / introduced in eq. ([2]) is expected to be smooth and universal, 
once the non-universal amplitudes ao and bo are fixed. 

2.2. Influence of the boundary conditions 

We now discuss those properties which do depend on the b.c. of the system. Extending 
the corresponding specific discussion in section [TJ here we shall consider the following 
b.c. 3 

• SP: semi-periodic b.c, i.e., periodic b.c. in both directions corresponding to the 
length L\\ and various b.c. in the direction corresponding to the length L, 

• FP: fully periodic b.c. 

In the case of SP b.c. two surfaces L» x Ln are present for which we shall consider 
either free or fixed b.cjj. Unless explicitly stated otherwise, the results presented in the 
following hold for both homogeneous and inhomogeneous b.c. 

| Here we shall not consider geometries with corners. Their presence gives rise to additional logarithmic 
singularities oc ln L in the free energy; see Ref . |68j and references therein for a discussion. 
§ Fixed boundary conditions can be obtained by applying infinitely strong surface fields. Here we do 
not address the dependence on scaling variables generated by the presence of finite surface magnetic 
fields, which results in interesting crossover phenomena |43U47| . 



F(t, H,L^oo, L,| -> oo) = /buik(t, #) + T / surf (t, H) + 0{e~ L /S/L), 
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As we mentioned in section 12. 1[ away from criticality the free energy density can be 
expanded into specific geometric contributions. For SP b.c. with laterally homogeneous 
(i.e., translationally invariant) b.c. on the surfaces and away from criticality, the free 
energy density J 7 has the following expansion as a function of the reduced temperature 
t and the external bulk field H [T61I22] : 

1 

L" 

where fbuik{t, H) is the bulk free energy density in the thermodynamic limit whereas 
/ sur f (t, H) is the surface free energy density associated with the surfaces of area L» x L\\ 
and L\\ — > oo. If the b.c. on the two surfaces are not translationally invariant in lateral 
directions, additional terms in the expansion of eq. (jSj) arise. In particular, in the 
presence of stripes with alternating b.c. an additional contribution oc 1/(LL\\) = p/L 2 
appears, which corresponds to a line free energy. In the case of FP b.c. / sur f = so 
that the free energy is simply proportional to the volume LL 2 , with exponentially small 
corrections. As anticipated in section I27TJ in eq. (JSJ) /bulk is independent of b.c. whereas 
/surf depends on the b.c. and on the local geometry, here taken to be planar. This allows 
one to identify a surface free energy for all types of surfaces, and the quantity f SU rf{t, H) 
in eq. (JSJ) is the sum of the two independent surface free energy densities associated with 
the confining walls of the corresponding semi-infinite systems [5],EB] ■ 

In the critical region the free energy is decomposed into singular and non-singular 
contributions (see eq. (0Q)). As in eq. (jHJ), the non-singular background jF( m > (t, H, L, p) 
can be decomposed further into geometric contributions 

7^{t, H, L -> oo, L,| -> oo) = f™{t, H) + I/£3(t, H) + 0( e - L /L) } (9) 

where the same considerations concerning the dependence on geometry and 
boundary conditions apply as explained after eq. (JHJ). In eq. (|HJ) the correction terms 
are characterized by a decay length of the order of the lattice constant, i.e., 1. In the 
following we shall neglect such corrections. 

Finally, also the universal scaling function /, which we have introduced in eq. (j2]), 
depends on the specific b.c. of the system. 

2. 3. Bulk, surface, and excess free energies 

Off criticality, the expansion in eq. (jHJ) applies in that form if the size of the system 
is large compared with the correlation length, i.e., in the limit L,Ln — > oo at fixed 
t ^ 0. This expansion defines the functions fbuik{t,H) and f BW t(t,H). Once this limit 
is taken, one can consider the limit \t\ — > [68] - On the other hand, off criticality 
the infinite- volume limit of the free-energy density can be calculated from the functions 
appearing in eqs. ([LIE]): in this limit the scaling variable tL yt diverges. In the absence 
of an external field H, the existence of a finite limiting value /buik(^ 5 H — 0) of the bulk 
free energy density requires that [72JI73] 

f(x, 0,p) ~ q±\x\ dv , \x\ -» oo, p < oo, (10) 
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with universal amplitudes q± and with the signs ± corresponding to t ^ so that the 
bulk free energy density decomposes into non-singular and singular terms: 

F{t -> 0, H = 0, L -► oo, L,| -)■ oo) = /£(£) + /^(t) = / Mk (t), (11) 

with the order of the limits given by L» — > oo, L — > oo, and £ — > 0. The singular part 
of the bulk free energy density is then given by 

/Sk(* ^0,H = 0)= J*\t, H = 0,L^oo,L ll ^ oo) 

= fSi^) = Q±\tr = q±\t\ 2 - a . (12) 

As we mention at the end of section I2.2[ the universal scaling function / depends 
on the b.c. However, the bulk free energy density is independent of the b.c. This 
implies that the asymptotic behaviour of / reported in eq. ffTUj) is independent of the 
b.c. Similarly, in the limit t — y the surface free energy density decomposes into 
non-singular and singular terms: 

f amt (t -y 0, H = 0) = f£>{t ^0,H = 0) + /<&(* -> 0, H = 0) (13) 

and the singular part of the surface free energies f^(t) is obtained by requiring that 
(compare eqs. (jHJ) and ([1])) the following limit exists: 

T T J™ < L (^ s \t, H = 0,L, L,|) - /«,(«)) = /«,(*). (14) 

This in turn requires that 

/ (x ->■ oo, p < oo) - g±|x|^ = r±|a;| , ' (d - 1) = r ± |x| 2 " as , (15) 

with «5 = a + z/ and universal amplitudes r±, so that the singular part of the free energy 
density becomes 

^ s \t,H = 0,L ^ oo,L|| -► oo) = f^(t) + y^(t), (16) 

with 

/£f(*) = r±|«r (d - 1) ) (17) 

which provides the leading correction to eq. (TT2|) . Since for FP b.c. there are no surfaces 
present, in eq. (fT5l) and eq. (TT7T) one has r± — and the singular part of the free energy 
density, in the large- volume limit, behaves as in eq. (fT2l) . For SP b.c. r± is the sum 
of the contributions from the two confining surfaces in the corresponding semi-infinite 
geometries. If one of these two surface exhibits a pair of individual chemical steps, 
its contribution is the mean value of those amplitudes which belong to those laterally 
homogeneous surfaces from which the chemical steps are made of. 

According to eq. (ITTj) the singular surface free energy density exhibits a singularity 
at the same temperature as the bulk free energy density. Concerning the critical 
behaviour at the surfaces, one has to distinguish different surface universality classes, 
which depend on the interactions at the surfaces. In particular, for the critical behaviour 
of an Ising system, the surface universality class depends on the so-called surface 
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enhancement cq, the sign of which indicates whether the order parameter at the surface 
is enhanced or diminished compared with its bulk value. The ordinary universality class 
holds for systems in which the interactions between the surface spins depress the order 
parameter and it corresponds to Co > 0. The extraordinary universality class holds if 
these surface interactions enhance the order parameter and it corresponds to c$ < 0. In 
this case the surface orders at a temperature higher than the bulk critical temperature. 
Still the bulk phase transition in the presence of an already ordered surface region 
results in a singularity of the surface free energy at T c . The case en = corresponds 
to the special universality class which belongs to a multicritical pointljjj. In the presence 
of surface fields, which explicitly break the Z2 symmetry, even in the absence of an 
additional bulk external field, the surface is always ordered. This corresponds to the 
so-called normal universality class which exhibits the same kind of singularities as the 
extraordinary universality class [THUS]. Finally, we note that fixed b.c. can be realized 
by infinitely strong surface fields, and thus they correspond to the normal/extraordinary 
surface universality class. We refer to Refs. [H,[5] for more detailed discussions of surface 
critical phenomena. 

The decomposition according to eq. (jSJ) and eq. (1T6|) provides a transparent 
interpretation only for ( <C L. In the critical region such a decomposition becomes 
blurred and the scaling function / in eq. (J5J) encodes all geometrical informations 
regarding the presence of surfaces, edges, etc. Under such circumstances it is possible to 
formally define a scaling function associated with the surfaces by comparing two system, 
with identical critical bulk behaviour but different b.c. [76]. To this end we consider a 
system with either FP or SP with homogeneous surfaces b.c; according to eq. (|HD, the 
free energy densities away from criticality are given by 

7- FP (t,if,L,L||) = /buik(i,#) + 0(e- L/ t/L), (18) 

and 

Jsp(£, H, L, L,|) = / bulk (t, H) + jf sw( (t, H) + 0(e- L ^/L), (19) 

where fbu\k(t,H) and f sm i(t,H) are independent of L and L\\ (compare eqs. (TT2l) and 
ffTTj) . respectively). We now introduce a new quantity 

/ surf (t, H, L, L\\)=L (F SP (t, H, L, L,|) - F FP (t, H, L, L„)) , (20) 

such that 

f suri (t,H,L^ oo,L|, -> 00) = f Bmi (t,H) + 0(e- L/ «). (21) 

By using eq. ([9]) and eq. © the singular part of / sur f(t, H, L, Ly) can be expressed as 

' / SP (u t L^ , u h W» , p) - / FP (utL* , u h L»> , p) 



&t,H,L,L ]{ ) L2 



(22) 



where / FP and / SP are the free energy scaling functions of J 7 ^ for FP and SP boundary 
conditions, respectively (see eq. (T5])). In contrast to the surface free energy density 

|| The special transition is located at Co = only within mean field theory; the fluctuations shift this 
point [S]. 
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/ S urf (t, H) in eqs. (J8]) and ( fT4|) . the extended surface free energy density / s ^ f (t, H, L, Ly) 
introduced in eq. ( 120}) depends not only on the local geometry but on the whole original 
reference geometry (only rinij, j £ || _ >00 / s ^ f reduces to f^). In particular for SP b.c. / s ^ f 
encodes the presence of both surfaces, the contributions of which cannot be disentangled, 
while as mentioned after eq. (jSJ), the surface free energy densities / s ^ f and / s ° r s f are the 
sum of the single- surface contributions in the corresponding semi-infinite geometries. In 
section 15.21 we shall extend this kind of reasoning in order to include line contributions 
to the free energy. 

Considering for simplicity a system in the absence of an external field H, the 
excess free energy /ix is defined as the remainder of the free energy density J 7 ^ after 
subtraction of the bulk contribution: 

f£(t,L,L ll )=^\t,H = 0,L,L ll )-fS lk (t). (23) 

According to eqs. (J2J) and (fT2l) it obeys the following scaling law: 

/£>(*, L,L|,) = ^A MAp). (24) 

The critical Casimir force Fc per area Lf, and per ksT is defined as 



1 c 



d(Lf, 



(8) 

ex 



dL 



t.L» 



(25) 



Analogous to eq. (J2J), in d = 3 the critical Casimir force exhibits the following scaling 
behaviour: 

F c (t,L,L Vi ) = jp0(u t ir,p) (26) 

with 

of \ i a -iw ^ dA(r,p) OA(t,p) 

0{r,p) = {d-l)A{T,p)-y t T — — p — , (27) 

so that at the critical point r = u t L yt = 

0(0, p) = (d- 1) A(0, p) - P^^- = O(p). (28) 

We now consider the Taylor series expansion at criticality 

oo 

A(0,p) = ^A nP ", (29) 

n=0 

which does not capture contributions which have a vanishing Taylor expansion, such as 
exp(— L/L\\) = exp(— 1/p). In this sense eq. (1281) yields a Taylor expansion for 0(p): 

oo oo 

©(p) = E( rf - l - n ^p n = E 9 *p n > ( 3 °) 

n=0 n=0 

such that 6 n= d-i =0. In d = 3 this implies that for p — y the critical Casimir force at 
the critical temperature has no quadratic term oc p 2 in its aspect-ratio dependence. 



Critical Casimir forces and adsorption profiles 10 

2-4- Corrections to scaling and improved models 

The scaling behaviour discussed in the preceding subsection is valid only up to 
corrections-to-scaling contributions. We distinguish two types of scaling corrections: 
non-analytic and analytic ones. The non-analytic corrections are due to the presence of 
irrelevant operators. In this case in eq. (j2]) an additional dependence on scaling fields 
arises, which are characterized by negative RG-dimensions. In the finite size scaling 
limit (FSS), i.e., for L — > oo, t — > at fixed £/L, this results in the following expression 
for the singular part of the free energy density J 7 ^ : 



^ s )(£,tf,L,Ln) = -i 



(31) 



/ (utL* , u h W» ,p)+J2 LkV% ^ (^tL yt , u h V>» , p) 

i,k>l 

where yi < 0, i > 1, are the RG-dimensions of the irrelevant operators and gt are smooth 
functions which are universal up to a normalization constant. The leading correction 
is given by the operator that has the smallest negative dimension. This is usually 
denoted with u, so that the leading scaling corrections are oc L~ u . For the standard 
three-dimensional Ising model one has u> = 0.832(6) |66j. 

In a family of models characterized by an irrelevant parameter A, it can occur that 
for a certain choice of A the leading correction-to-scaling term oc L~ w vanishes: such a 
model is called improved. In such models, the observed scaling corrections usually decay 
much more rapidly: as L~^ 2 with u>2 ~ 1.67(11) [77] for the three-dimensional Ising 
universality class. Improved models have turned out to be instrumental for obtaining 
high-precision results for critical phenomena [67]. In a lattice model there are also 
scaling corrections oc L~ WNH due to the presence of non rotationally-invariant irrelevant 
operators, i.e., operators which break the rotational invariance. The leading one for 
lattices with cubic symmetry leads to cjnr — 2 [78] . 

Another type of scaling corrections is provided by the so-called analytic scaling 
corrections, which can stem from various sources. The linear expansion of the scaling 
fields introduced in eqs. ([3]) and (j4]) is valid only up to higher-order terms in the 
expansion. For instance, it is easy to see that an additional term oc t 2 in eq. ([3]) gives rise, 
in the FSS limit, to scaling corrections oc L~ l l u . It is worthwhile to note that, for the 
Ising universality class and general O(N) models, one has v < 1 so that these corrections 
are usually not observed because they are subdominant. Analytic corrections can be 
also due to boundary conditions: not fully periodic b.c. induce additional corrections, 
which are proportional to L~ 1 , It was first proposed in Ref. [79j in the context of 
studying surface susceptibilities, that such scaling corrections can be absorbed by the 
substitution L — > L + c, where c is a non-universal, temperature-independent length. 
This property has been recently checked numerically in Refs. (BTJ [HTil El] for the XY 
model for SP b.c. with free surfaces and in Ref. [M] for the Ising model for SP b.c. with 
homogeneously fixed surface spins. We note that, in the same spirit of the expansion 
of eqs. ([3]) and (HI), one should expect this substitution to be correct up to higher 
orders in 1/L, i.e., L — > L + c + 0(1/L). However, such higher-order terms would 
induce corrections to scaling oc L~ 2 which interfer with those due to the aforementioned 
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breaking of rotational invariance; accordingly, scaling corrections oc L~ 2 are captured 
by the more general ansatz given by eq. fl3TT) . 

In the present contribution we study the critical Casimir force using an improved 
model in three spatial dimensions, employing SP b.c. with fixed surface spins. On 
the basis of the above discussion for such a model the leading scaling corrections are 
expected to be proportional to L~ x . Furthermore, we assume that also in this case in 
leading order we can absorb such a scaling correction by the substitution L — > L + c. 
In this sense eq. ( 126]) is replaced by 



F c (t, L, L,|) = jj^—^0 (a t(L + c)« ±±^ , (32) 

where we have used the expansion of eq. (E]) to leading order, thus neglecting corrections 
oc L~ x l v ~ L~ L6 . We have also neglected corrections to scaling due to next-to-leading 
irrelevant operators, which are expected to scale as L~ U2 , U2 ~ 1.67(11) [77]. If we 
expand eq. (132]) for L — > oo, at fixed tL Vt and p, we obtain the following expression for 
the leading behaviour and for the scaling corrections of the critical Casimir force: 

F c (t,L,L\\) = -1 [e{a Q tLy\p) + li/>(aotL*,p) 

if \ - oaf \ , l d9 ( T 'P) , d °( T iP) /oon 

$\T,p) = -39[T,p) + -T — ^ + P — ^ • 33 

V OT Op 

The comparison of eq. ( 13"3"|) with eq. ( 13T|) reveals an important difference in the 
correction-to-scaling terms. While in eq. ( 13T]) the functions gi, which provide the 
scaling corrections, are independent from the leading scaling function /, in eq. ( 1331) the 
correction-to-scaling function ip follows from the leading scaling function 9. In principle 
this prediction provides an opportunity to check the appealing ansatz in eq. ( !32l . 

2. 5. Finite-size scaling of the order parameter 

The scaling properties of the order parameter can be determined from a generalization 
of eq. (EJ) by introducing a spatially varying external field (HE]. Here we consider SP b.c. 
with fixed surface spins. The system is described by the coordinates x, y, z on the lattice, 
with z the coordinate corresponding to the length L. Furthermore we restrict ourselves 
to the case in which the system is translationally invariant in the y direction, i.e., the 
b.c. of fixed surface spins are translationally invariant in one direction, corresponding 
to the length L\\. The order parameter $ (i.e., the magnetization per volume) exhibits, 
in leading order, the following scaling form: 

$(t,x,z,L,L ll ) = B\tfP ± (^,-^,^,p), £± = ^\tr, (34) 

where £± is the correlation length for t ^ with non-universal amplitudes £ (see 
eq. flSJ)). The scaling functions P±(x,z,L, p) (for t ^ 0) are universal and depend 
also on the universality classes of the two confining surfaces. B is the non-universal 
amplitude of the spontaneous bulk magnetization at T < T c : 

$buik(t -► 0") = B\tf. (35) 
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The amplitude of the scaling functions P± is fixed by the requirement that eq. (1551) 
is reproduced in the bulk limit, i.e., hm.£_ >oo P_(:r, z = L/2,L,p) = 1. Due to the 
fixed b.c, in the limits t — y O^ 1 the scaling functions P± give a unique, nonvanishing 
order parameter profile, respecting the universality of the ratio £tT/£o~ — 1-896(10) [7T] . 
Equation (1341) can be rewritten as 

/ L \ _/3/l/ f x z L \ 
${t,x,z,L,L ll )=B( T ±) ± (-,-,_, p ), (36) 



sv "vi'i'e 



;± 



where 



0±(x, z, Z, p) = L^Pi (x = xL, ~z = zl, L, p) (37) 

are also universal scaling functions. 

If the order parameter on the whole of each surface is either fixed to +1 ((+) b.c. 
in the following) or —1 ((— ) b.c.) the system is translationally invariant in both lateral 
directions. In this case the order parameter does not depend on x and eqs. (1511) and 
(136]) reduce to 



^t,z,L,L {l ) = B\tfPt^(^-,^p ) j 



T \ ~M V / , T 

i-' \ , (horn"! I Z 1j 



(38) 

<&(*, Zl L, L,|) = B (^ 1 0i hom) (|, ^, p) , (39) 

with, again, universal scaling functions P± (2, L,p) and <p± (z,L, p) = 
LP/ u P± m (z = 2L,L,p); the superscript (horn) indicates laterally homogeneous b.c. 

For fixed t, fixed z, fixed aspect ratio p (p < oo), and for increasing separation L 
the scaling functions introduced so far converge to the corresponding half-space (semi- 
infinite) scaling functions for the confining surface located at z = 0: 

P±(x,z,L^ oo, p) ->P± l0 o(x,z), (40) 

pf om) (5,L^oo,p)->p£ m) (5), (41) 

where P± )00 (x, z) and P± &(z) describe the order parameter profiles for a semi-infinite 
(oo/2) geometry in the presence of an inhomogeneous and a homogeneous surface, 
respectively, at z = 0: 

${t,z) = B\t\ /} P£g? (-?-) horn, surf., oo/2, (42) 



$(£, x, z) = B\tfP ±iOC ( -^ , — ) inhom. surf., oo/2. (43) 



X z 

In the limit L — y oo at a fixed aspect ratio p, one also has L» —y oo. Accordingly, 
the limits in eqs. (14*01) and (141j) describe a system with an infinite lateral extension. 
In such a geometry no quantity depends on the original aspect ratio p of the confined 
system. Therefore we have dropped the p-dependence on the r.h.s. of eqs. ( 1401) and ( 14T1) . 
as well as in eqs. ( 1421) and ( 14"3]) . For homogeneous surfaces in the film geometry, the 
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short-distance behaviour of the scaling functions in eq. (I3"8j) and (J121 is given by [5]I8"2"] 

p |hom) ( - ^ 0> ^ ^ _ C± ^-^, (44) 

Pi h ™\z^0)^c ± z~^, (45) 

i.e., the leading behaviour near the surface is not influenced by the finite film thickness 
which, however, gives rise to distant-wall corrections (see below )% Equations (l4~4"l) and 
( 145]) imply that in both cases at criticality the order parameter varies algebraicly: 

$(t = 0, z -> 0, L, L,|) ~ 5c±(^ ± )-^. (46) 

The amplitudes c + and c_ are universal [H2] with c+/c = (£o~/£o~) -/3 ^- The presence 
of a second wall induces corrections [83,84j to the leading decay (eq. (14"B"]) ) of the critical 
order parameter. For Lii — > oo the leading distant-wall corrections are given by 

<f> ab (z, t = 0) = $ a>00 (z, £ = ti)?ML) (47) 

with 

J-«6(z/L) = 1 - C a (d- l)Aa6 (£)* , (48) 

where $ aj b is the order-parameter profile in the presence of walls with surface universality 
classes a and h and <& a ,oo is the semi-infinite order-parameter profile for the near wall, 
given by eq. (142]) . In the universal scaling function jF^z/L) the coefficient C a depends 
only on the boundary conditions of the wall a, whereas the influence of the boundary 
condition at the distant wall enters via A ab = A a b(0, 0) (eq. ff24|) ). which is the Casimir 
amplitude of the free energy at criticality. 

The expressions in eqs. (EI}, (ESJ, (EH}, flU}, ©, and g3J are valid up to 
corrections to scaling. In analogy to eqs. ( 13T|) and (|32|) we expect, to leading order, the 
occurrence of corrections to scaling oc L~ w due to irrelevant operators and of corrections 
to scaling oc L^ 1 due to the boundary conditions. Since we shall study an improved 
model, the latter ones are expected to be the leading corrections to scaling. 

3. Model and method 



As discussed in section 12.41 in order to have a better control on the finite-size scaling 
corrections, we resort to an improved model belonging to the 3D Ising universality class. 
To this end we chose the Blume-Capel model |85U86j . It is defined on a three-dimensional 
simple cubic lattice, with a spin variable Si on each site i which can take the values 
Si = — 1, 0, 1. The reduced Hamiltonian for nearest neighbour interactions is 

■H = -(3j2 S > S J + D Y. S i> St = -1,0,1, (49) 

(*3) i 



% In principle, the r.h.s. of eq. (T44|) could depend on the aspect ratio p. However, since a confined 
system displays a finite correlation length and the order parameter is a local quantity, we expect that 
the effect of a nonzero aspect ratio p, i.e., a finite lateral size, will be effectively very weak, if not absent 
at all. Thus for simplicity in eq. (J44I) we have omitted a possible dependence of the constants c± on p. 
The data which will be presented in section [5] support this observation. 
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so that the Gibbs weight is exp(— %). In line with the convention used in Refs. [64"jl7ip87j . 
in the following we shall keep D constant, considering it as a part of the integration 
measure over {Si}, while we vary the coupling parameter /3, which is proportional to the 
inverse temperature, /3 ~ 1/T. The relation between the dimensionless parameters (3 and 
D characterizing the reduced Hamiltonian in eq. ( 1491) and the physical parameters can be 
obtained by introducing coupling constants J and D, so that the reduced Hamiltonian 
at temperature T reads 



n = — ( -jJ2 s ^ +f) J2 s * I > ( 5 °) 

(ij) i 



k B T 



where J, D, and k B T have the dimension of energy. The comparison of eqs. (|49p and 
( 150]) gives (3 = J /(ksT) and D = D/(ksT). In the limit D — > — oo, one recovers the 
usual Ising model, because in this limit any state for which there is an io such that 
Si = is suppressed relative to the states {Si = ±1}. For d > 2, the model exhibits a 
phase transition at (3 C = (3 C (D) which is second order for D < D tTi , and first order for 
D > D tri . The value of D tri in d — 3 has been determined as D tli ~ 2.006 [88J and as 
D tri ~ 2.05 [89J. At D = 0.641(8) [87J the model is improved, i.e., leading corrections 
to scaling oc L~ w with u = 0.832(6) [66JJJ are suppressed. At this value of the reduced 
coupling D the model is critical for — (3 C — 0.3856717(10) [71]. We mention that for 
this improved model, i.e., at D = 0.641, 25 th -order high-temperature expansion series 



are available [7T]. In Appendix A from these series we infer (in units of the lattice 
constant) 

£+ = 0.415(2) (51) 

as the value of the non-universal amplitude of the true correlation length above T c (see 
eq. ©). 

In the following we use as the values of the critical exponents v = 0.63012(16) and 
7] = 0.03639(15) which have been obtained by analyzing the 25 th -order high-temperature 
expansion series for three improved models within the Ising universality class [71] . 

In order to determine the critical Casimir force, we follow the approach introduced 
in Ref. [59J, which we briefly describe here. For two reduced Hamiltonian ~Hi and I-L2 
associated with the same configuration space {C} we construct the convex combination 
H(X) 

n(X) = (l-X)H 1 + XH 2 , AG [0,1]. (52) 

This Hamiltonian H(X) leads to the free energy F(A) in units of k B T. Its derivative is 

dF(A) 1^ {C] ^x~ e r „N 

ox j: {c} e- m) ' [ } 

+ This paper provides an updated value of the coupling D — 0.656(20) for which the model (eq. ([49]) ) 
is "improved" [57] . as well as an updated value of /3 = 0.38567122(5) for which the model is critical at 
D = 0.641. 
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Figure 1. Film geometry with aspect ratio p = L/Lu confined by a homogeneous 
upper surface and by a lower surface with a chemical step. 



Combining eqs. ( 152]) and ( )53|) we can determine the free energy difference as 



F(l) - F(0) 






d\(U 2 - Hi); 



(54) 



where (H2 — T-i\)\ is the thermal average of the observable H2 — Hi with the statistical 
weight exp(— H(A)). For every A this average is accessible to standard Monte Carlo 
simulations. Finally the integral appearing in eq. ( l54"j) is performed numerically, yielding 
the free energy difference between the systems governed by the Hamiltonian % 2 and Hi, 
respectively 

To be specific, we consider a three-dimensional lattice Lx L\\x L\\ with SP b.c. (see 
section I27L]) . with periodic b.c. in the lateral directions x and y, and fixed spins at the 
two surfaces z = and z — L — 1, so that there are L — 2 layers of fluctuating spins. 
The spins at the upper surface z = L — 1 are fixed to +1, and on the lower surface z = 
we employ a single, straight chemical step, where the surface is divided into two halves, 
one (x < 0) with spins fixed to —1 and the other half (x > 0) with spins fixed to +1. 
The presence of lateral periodic b.c. in the x direction generates an additional chemical 
step at the lateral boundaries, resulting in a system with a pair of individual chemical 
steps. This geometry is illustrated in figure [TJ 

Following Refs. (59j[60], we apply eq. ( 154")) with Tii as the Hamiltonian of the lattice 
shown in figure [1] and I-L2 as the Hamiltonian of a lattice (L — 1) x L11 x L\\ plus a 
separate two-dimensional layer, so that both Hamiltonian share the same configuration 
space. This layer can be inserted into the film by varying the coupling (1 — A)/3 with its 
neighbouring planes between and /3. With this we evaluate the following quantity: 

1 ^ 

a 



1 r 1 

J(/3,L,L||) = -2 / d\(H 2 -Hi); 

L \\ Jo 



(55) 



In accordance with eqs. (jSJ) and ( )23|) . the free energies (in units of ksT) corresponding 
to %i and %2 can be expressed as 



Fi = F(0) 



LLlf, 



(ns), 
bulk I 



-2 /(ns). 



? (ns) 



2Ab) 



(*) + HfSi® + hfZ'it) + LL lf Z k (t) + LLifS(t, L, L„), (56) 
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F 2 = F(l) = {L- l)L\f^l{t) + LJj/2(«) + L||/£?(t) 

+(£ - l)4/bl(*) + (L ~ 1)4/W(*, ^ - l^ll) + 4/^W, (57) 

where /^ (£) represents a (possible) line contribution to the non-singular part of the 
free energy due to the presence of the pair of individual chemical steps (see the discussion 
after eq. (jHJ)) and /2d(£) is the free energy density per area L 2 and in units of ksT of 
the additional 2D layer, which is not critical near the phase transition in d = 3 we are 
interested in. By combining eq. (155]) with eqs. (1541 . ( 1561 . ( |5Tj) . and (ITT]) one has 

/ (/?, L, L,|) = [f 2D (t) - f hulk (t)} - [/«(*, L, L,|) - /*?(*, ^ - 1. £||)] - (58) 
with 

w«) = /££(*) +/&(«), (59) 

fi>(i,L,L|l)=L/M(t,L,L|l), (60) 

*-^- (61) 

A Taylor expansion of the last term in eq. (158]) around the film thickness L — 1/2 gives 

mt,LM) - f£(t,L - 1,L„) = ^f (*, L - i L„) + ^^ (t, L - i,L„).(62) 
By combining eqs. (160]) and (1251 . we have 

^(t,L,L||)=-F c (t,L,L||), (63) 

i.e., the first term in the expansion (162]) corresponds to the critical Casimir force. 
Neglecting for the time being corrections to scaling, we substitute eq. (125]) in eq. (JHSD 
obtaining 

-^L(t,L,L ll ) = -—6(u t L«*,p) = -—9(T,p), T = u t V% y t = -. (64) 

The higher order derivatives follow from eq. (I64|) : 

d 2 f ^ 3 v 80 1 89 1 

-qJT- (*, L > L ll) = 71^ ( r ' P) ~ JZ T fr ( r > ^ " -rJP-Q-p (T ' P) " T^ 2 (r ' P) ' (65) 

-^T (*> L ' L ll) = "J/* (r, p) + -^r— (r, p) + — p— (r, p) = —0 3 (r, p) . (66) 
Inserting eqs. (1631 . (164]) and (166]) in the expansion (1621 we obtain: 
/£)(*, L,L,|) -/£)(i, L-1,L|,) 



-A', (t,L-i L| 



1 3 («t(i-l/2) w ,p) 



24(L-l/2) 2 ^K(L-l/2r,p) • (6?) 

In the FSS limit, i.e., in the limit L — > 00 at a fixed ratio £/L or fixed r, the term in 
brackets represents a scaling correction ~ L~ 2 . Such a correction is negligible relative 
to the leading scaling correction ~ L~ x which we expect for the model under study here 
(see the discussion in section 12741 . It is easy to see that higher order terms in the Taylor 
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Figure 2. Decomposition of a system with a chemical step at the lower surface in the 
limit of vanishing aspect ratio p = L/Lii — > 0. The symbol + (— ) indicates regions in 
which the boundary spins are fixed to +1 (—1) (see the main text). Periodic b.c. are 
imposed on the lateral directions. 



expansion of eq. (162]) result in additional corrections ~ L , k > 1, which are negligible 
as well. 

By inserting eq. fl67|) in eq. f )58|) and neglecting corrections ~ L~ 2 we finally obtain 



/(/3,L,L,|) = [/ 2D (i)-/ bulk (t)] + F ( 



c 



t,L--,L\ 



(68 



which still involves the subtraction of the bulk free energy density and the areal free 
energy density of the two-dimensional layer. Upon substituting eq. (132|) into the previous 
expression, we expect the following scaling form for the quantity in eq. (l68l) : 

9 (a Q t(L - 1/2 + a) 1 /", (L - 1/2 + c)/L\\) 



I{P,L,L\\) =B(t) 



B(t) 



(L - 1/2 + c) 



9(t((L-l/2 + c)/ti) 1/u ,p(l + (c-l/2)/LJ) 



(L-l/2 + c)- (69) 

where B{t) = f2Dif) — /buik(^) is a L-independent background term, and 9(r, p) is the 
universal scaling function associated with the critical Casimir force. It is important to 
notice that simulations of a lattice with aspect ratio p = L/L\\ provide data for the 
scaling variable corresponding to a different aspect ratio 

i/r 



p(L) =p[l 



(70) 



which of course converges to p for L — )■ oo. 

In the following, we consider the critical Casimir force for various aspect ratios 
p, in particular also the extrapolation to the slab limit, i.e., for p — > 0. In this limit 
the critical Casimir force reduces to the mean value of the critical Casimir force in the 
++ configuration, for which all boundary spins are fixed to the same value, and of the 
critical Casimir force in the H — configuration, for which at one surface all spins are 
fixed to +1 and at the other surface all are fixed to —1. This can be understood by the 
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following argument. Near the bulk critical point T c the actual correlation length in a 
slab is bounded by L (it may diverge laterally at the film critical point T c j(L) < T c ). 
Therefore near T c , in the presence of a pair of individual chemical steps, one can divide 
the system into four pieces: two pieces around each chemical step, which are influenced 
by the steps, and two lateral pieces which, sufficiently far from the steps, are in the ++ 
and H — configuration, respectively (see figure [2] for an illustration). The former two 
pieces have a lateral size ~ L, while the latter ones have a size ~ L»/2 — L. Thus the 
total critical Casimir force per area is given by 

F c ~ -^ [F ++ L||(L||/2 -L) + F+_L||(L||/2 - L) + 2F p L ll L] , (71) 

ii 

where, in accordance with figure |2j F ++ is the critical Casimir force per area in the 
++ configuration, F + _ is the corresponding one in the H — configuration, and F p is the 
force per area due to central and lateral pieces. By taking the limit L\\ — > oo at fixed 
L, or equivalently, the slab limit p — y at fixed L, one finds Fc = (F ++ + F + _)/2 Q 
Along the same line of reasoning one can consider a system in which both confining walls 
exhibit a chemical step. In the case that equal (opposite) b.c. face each other, for p — > 
the critical Casimir force is expected to be the same as for laterally homogeneous ++ 
(H — ) b.c. on both sides. This holds even if the two chemical steps are shifted relative 
to each other. In the case of equal b.c. facing each other and in the presence of a bulk 
magnetic field the force again equals the mean value of the ones for the two halves. 

4. Critical Casimir amplitude at T c 

As discussed in the preceding section, the quantity l(/3,L, L\\) can be computed by 
standard Monte Carlo simulation techniques combined with a numerical integration. 
At the bulk critical temperature eq. (169]) becomes 

/(ft,L, i|l )=B + ™, (72) 

{L + c) 6 

where B = B(t = 0) is the non-universal background, O (p) = 9 (0, p) is the universal 

amplitude of the critical Casimir force at T c , L = L — 1/2, and p(L) is given by eq. ( 17D1) . 

In a series of Monte Carlo simulations, we have evaluated the quantity / (f3 c , L, Lu) 

for lattice sizes L = 8, 10, 12, 16, 20, 24, 32, 40, 48, 64 and aspect ratios p = 1/6, 



1/8, 1/10, 1/12. Certain details of the simulations are reported in Appendix B From 
eq. (J72l we can obtain an estimator for the critical amplitude 0(p) by considering the 
difference concerning two different lattices. To this end we introduce the quantity 



, 3 /(&,L,L|l)-j(&,aL,aL|l) 

l-((L-l/2)/(aL-l/2)) 3 ' 



e est (L, P ) = {L- 1/2) 3 : yc \,: 'Jr: zi ^) 



with a fixed integer a. Using eq. (172")) . one finds that Q es t(L — >■ oo,p) — > 0(p), with 
corrections ex 1/L. Moreover such corrections are proportional to c and to p, i.e., in the 



* 



We thank Ettore Vicari for pointing out this argument. 
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Figure 3. The estimator Q cst (L,p) (eq. (|73)0 of the critical Casimir amplitude at 
T c , for aspect ratios p — 1/6, 1/8, 1/10, 1/12. The simulation data correspond to 
the system shown in figure [T] for the improved Blume-Capel model (|49j) and for the 
standard Ising model (Is.). The statistical error bars of the Monte Carlo data are as 
large as the symbol size or smaller. 



slab limit p — > and in absence of the leading scaling corrections (c = 0) the estimator 
given in eq. (I73"j) equals 6(0) up to corrections oc L~ W2 , uj 2 = 1.67(11) [77]. We show the 
quantity Q est (L,p) in figure [31 for a = 2, as obtained by a simulation of the improved 
Blume-Capel model ( H9|) with D = 0.641. As a comparison, we also display some data 
obtained for the standard Ising model corresponding to the Blume-Capel model in the 
limit D —7- — oo. 

The prediction that for the improved model the leading corrections are proportional 
to 1/L is confirmed by the data shown in figure [31 Moreover the scaling corrections 
appear to be numerically small: the data at L = 8 differ from those at L = 32 by 3%, 
and the amplitudes for L > 20, i.e., 1/L < 0.0513 are all compatible within error bars. 
On the other hand, the data for the Ising model exhibit stronger corrections; for this 
model, we expect additional corrections oc L _aJ , with u = 0.832(6) [66] . 

In order to obtain an estimate of the critical Casimir amplitude 0(p), we expand 
eq. ( T72l) in terms of 1/L at a fixed aspect ratio p. In lowest order in 1/L one obtains 

©(p) , C(p,c) 



J(/3 C ,L,L||) =B 



I 3 



I 4 



L = L-l/2, 



(74) 
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Table 1. Fit of Monte Carlo data to eq. (|74|) with free parameters B, 9, and C. 
£min is the smallest lattice size taken into account for the fit. DOF denotes degrees 
of freedom. 



p=l/6 



p=\/l 



p = 1/10 



p = l/l2 



X 2 /DOF = 6.5/7 
B = 0.03235167(9) 
6 = 2.048(3) 
C = 0.60(2) 



X 2 /DOF = 14.1/7 
B = 0.03235196(8) 
6 = 2.132(2) 
C = 0.55(2) 



X 2 /DOF = 12/7 
B = 0.03235157(10) 
6 = 2.190(3) 
C = 0.47(3) 



X 2 /DOF = 3.9/6 
B = 0.0323514(1) 
= 2.229(3) 
C = 0.41(2) 



10 x 2 /DOF = 6.1/6 
B = 0.0323517(1) 
6 = 2.046(4) 
C = 0.63(4) 



X 2 /DOF = 7.6/6 
B = 0.03235210(9) 
6 = 2.124(3) 
C = 0.65(3) 



X 2 /DOF = 7.3/6 

5 = 0.0323517(1) 

6 = 2.181(4) 
C = 0.57(4) 



X 2 /DOF = 2.9/5 
B = 0.0323515(2) 
= 2.225(4) 
C = 0.46(4) 
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X 2 /DOF = 5.3/5 
B = 0.0323517(1) 
9 = 2.050(6) 
C = 0.56(7) 



X 2 /DOF = 7.2/5 
B = 0.03235206(10) 
9 = 2.126(4) 
C = 0.61(5) 



X 2 /DOF = 6.9/5 
5 = 0.0323517(1) 
9 = 2.184(5) 
C = 0.52(7) 



X 2 /DOF = 2.5/4 
B = 0.0323516(2) 
9 = 2.221(6) 
C = 0.51(7) 
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X 2 /DOF = 5.3/4 
B = 0.0323517(2) 
6 = 2.05(1) 
C = 0.5(2) 



X 2 /DOF = 4.9/4 
B = 0.0323522(1) 
9 = 2.112(8) 
C = 0.9(2) 



X 2 /DOF = 6.6/4 
B = 0.0323517(2) 
6 = 2.18(1) 
C = 0.6(2) 



X 2 /DOF = 0.6/3 
B = 0.0323518(2) 
= 2.20(1) 
C = 0.8(2) 



where C(p,c) = — 3c9(p) + (c — l/2)p0'(p) is the leading correction-to-scaling term, 
which depends on p and c. (See also the discussion below leading to eq. (!86|) .) However 
eq. (174]) provides an unbiased determination of the critical Casimir amplitude 6 (p), 
which is correct even if the corrections to scaling oc 1/L do not have an analytic origin. 
We directly fit our Monte Carlo (MC) data for the quantity I (/3 C = 0.3856717, L, Ly) to 
eq. f !74|) leaving B, 6, and C as free parameters. In order to control a possible systematic 
error due to subleading scaling corrections, we repeat the fit disregarding the smallest 
lattices. For the various aspect ratios in table [1] we report the fit results as a function 
of the smallest lattice size L min taken into account for the fit. 

Inspecting the fit results, we generally observe a good x 2 1 'DOF {DOF is the number 
of degrees of freedom, i.e., the number of statistically independent points minus the 
number of fit parameters) for L min > 10 and the results appear to be stable with respect 
to L min . While there is a clear dependence of the Casimir amplitude 6 on p, as expected 
(see eq. (169]) ) the background term B does not exhibit a dependence on p. We observe 
that B as determined at p = 1/8 is slightly shifted with respect to the corresponding 
values determined for the other aspect ratios. However the difference is tiny (two error 
bars in the worst case). Thus we conclude that the observed shift can be interpreted as 
a statistical fluctuation; note also that the ratio x 2 1 'DOF is slightly worse for the data 
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Figure 4. Critical Casimir force amplitude 6(p) = 6(0, p) (see eqs. (|26j) and (|28|) ) at 
T c for aspect ratios p = 1/6, 1/8, 1/10, 1/12, as inferred from tabled] (eqs. ([75 ]) -([75 ]l ). 
as well as the limit p — !> reported in eq. (|82[) . The dashed line represents eq. ()79)) . 
with its parameters given by eqs. (|82|) and (f83|) . As expected on the basis of eq. 130|) . 
0(p — >■ 0) does not exhibit a quadratic term. The statistical error bars have the same 
size as the symbols. 



at p = 1/8. 

Corrections to eq. ( |74l) are generated by next-to-leading irrelevant operators and 
result in an additional term oc L~ 3 ~ W2 , with uj2 = 1.67(11) [77]. Fits including such 
correction do not result in significant deviations from the results given in table [TJ 

By judging conservatively the variation of the resulting 9 with respect to L min , 
from table [1] we obtain the following estimates: 



0(1/6) =2.048(6), 
6(1/8) =2.126(5), 
6(1/10) = 2.183(6), 
6(1/12) = 2.223(7). 



(75) 
(76) 
(77) 
(78) 



These amplitudes are shown in figure HI They clearly show a linear dependence on 
the aspect ratio p. This allows us to determine the limit 6(p — > 0) by expanding 

6(p) = 6(0) + Ep + 0(p 3 ), (79) 

neglecting terms which are nonlinear in p. As discussed at the end of section 12.31 in the 
expansion of eq. (1791) . the amplitude of the critical Casimir force at T c has no quadratic 
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Table 2. Fits of Monte Carlo data to eq. ff80]) . with free parameters B, 6(0), E and 
setting c = 0. L m - m is the smallest lattice size taken into account for the fit. DOF 
denotes the degrees of freedom. 

6(0) E X 2 /DOF 



8 0.03235066(4) 2.4366(9) -2.026(8) 883/36 

10 0.03235100(4) 2.431(1) -2.05(1) 342/32 

12 0.03235128(4) 2.425(2) -2.07(2) 133/28 

16 0.03235153(5) 2.415(3) -2.06(3) 59/24 



Table 3. Fits of the Monte Carlo data to eq. ([80]) with free parameters B, 6(0), E, 
and c. i m in is the smallest lattice size taken into account in the fit. DOF denotes the 
degrees of freedom. 



6(0) E c X 2 /DOF 



8 


0.03235170(5) 


2.386(2) 


-1.991(8) 


-0.058(2) 


56/35 


10 


0.03235181(5) 


2.383(2) 


-2.02(1) 


-0.069(3) 


35/31 


12 


0.03235180(6) 


2.388(3) 


-2.05(2) 


-0.066(5) 


32/27 


16 


0.03235188(7) 


2.377(6) 


-2.04(3) 


-0.09(1) 


27/23 



term oc p 2 . Thus corrections to the above expression are expected to be at least of third 
order, i.e., oc p 3 . By inserting eq. (!79|) into eq. ( 1721) and using eq. (ITU]) we obtain in 
linear order in p 

I ^ L ^ ) = B + m±MLt^BI^iim, TsL _ 1/2 . (80) 

(L + c) 
Inspection of eq. (l80l) shows that scaling corrections oc 1/L for I emerge from two 
contributions. The first contribution stems from the dependence on the aspect ratio p: 
the MC data of a system with sizes L x L\\ X L\\, corresponding to p = L/L», result in 
the force for a system with perpendicular size L = L — 1/2, which has an aspect ratio 
(L — l/2)/L|| = p(l — 1/(21/)). This correction is proportional to p. The second source 
of scaling corrections is due to the non-periodic boundaries and it is proportional to c. 

We first consider fits of the MC data ignoring scaling corrections due to the 
boundary conditions. We fit all data for the various aspect ratios to eq. ( l80l) . setting 
c = and leaving B, O(0), and E as free parameters. The results for those fits are 
reported in table [2j 

These fits have a large ratio x 2 /DOF and in addition the fitted parameters show a 
systematic drift, which is larger than the statistical error. From this we conclude that 
scaling corrections due to the parameter c in eq. (IHOl) are sizeable within the statistical 
precision of the MC data. In view of this we fit all data for the various aspect ratios 
to eq. (|80|) . leaving B, 6(0), E, and c as free parameters. The results of these fits are 
reported in table El 

Here the situation is much improved and fits with L > 10 show a good ratio 
X 2 /DOF. Moreover, the results are stable upon increasing the smallest lattice size 
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taken into account, which underscores the quality of the fits. Accordingly we obtain as 
final estimates 

5 = 0.03235181(7), (81) 

6(0) = 2.386(5), (82) 

E = -2.04(3), (83) 

c= -0.066(5). (84) 

Also in this case we have performed fits of our Monte Carlo data by adding a next- 
to-leading scaling correction term to eq. flHOj) . These results do not exhibit significant 
deviations from eqs. f ]8~T]) -( j8"4"]) . In figure|4]we compare the amplitude of the Casimir force 
at T c for the various aspect ratios reported in eqs. (175l-(178]) with the linear dependence 
on p predicted by eq. ( 179]) . using eqs. ( 182]) and (183]) . We find very good agreement. 

As discussed at the end of section |3j in the limit p — > the critical Casimir force 
reduces to the mean value of the forces valid for the ++ and the H — configuration, 
respectively, where the first case corresponds to a system in which all boundary spins 
are fixed to the same value, and the latter case corresponds to a system in which 
on one surface the spins are fixed to +1 and on the other to — 1. These laterally 
homogeneous Ising systems have been investigated in Refs. [59]i6Q]|6l] by Monte Carlo 
simulations. According to Ref. [60], ++ = —0.76(6) and + _ = 5.42(4), so that 
(0 ++ + + _)/2 = 2.33(4), in marginal agreement with eq. ( 182]) . According to Ref. [64], 
++ = -0.820(15) and 6+_ = 5.613(20), so that (0 ++ + 6 + _)/2 = 2.396(13), in 
perfect agreement with eq. (1821 . In Refs. [59],[60] simulations have been carried out for 
the standard Ising model, for which the scaling corrections proportional to L~ 08 and 
those proportional to L^ 1 are difficult to disentangle. Thus we expect our results to 
be more reliable with respect to those of Refs. [59j|60], because in our model scaling 
corrections are under control. Moreover, our results are in perfect agreement with those 
of Ref. [64] , in which the improved Blume-Capel model has been used. 

Expansion of eq. (1801 to the lowest order in the scaling corrections yields, up to 
terms linear in the aspect ratio p, 

I (ft, L, h )=B + 9 <°L| Ef ' + - 3C9(0) ~ 2 S E - PE/2 ■ (85) 

L L 

Comparing this result with eq. (1741) and taking into account eq. (179]) . we obtain an 

expression for the amplitude C(p,c) which appears in eq. (1741 : 

C (p, c) = -3c6(0) - 2pcE - -pE. (86) 

Thus in line with the expression for C (p, c) given after eq. (174]) . the amplitude C (p, c) 
of the correction to scaling in eq. ( 174"]) can be expressed in terms of c and the critical 
Casimir amplitude O(p). This relationship is due to the analytic origin of the corrections 
to scaling (see also the discussion at the end of section I2.4p . Using the results in 
eqs. (Ig2> (154l) we obtain C(l/6) ~ 0.60(3), C(l/8) ~ 0.57(4), C(l/10) ~ 0.55(4), and 
C(l/12) ~ 0.53(4). Comparing these values with those shown in tabled] for the most 
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reliable data with L min = 10 and L min = 12 we observe agreement within the error 
bars. Given the limited available precision, this is a non-trivial consistency check of our 
scaling ansatz. 

5. Casimir scaling function 

5.1. General results 

In order to determine the full scaling function 9(t, p), in eq. ( 1691) we have to subtract the 
L— independent term B(t) from the quantity I(/3,L,L\\) sampled by MC simulations. 
Moreover, eq. (169|) tells that the MC data obtained at an aspect ratio p refers to the 
Casimir scaling function at a modified aspect ratio p given by eq. ( 170|) . Since we have 
determined the temperature-independent length c at T c (see eq. (184")) ). we can calculate 
the quantity p for all MC data. In order to subtract the background term B(t) we 
proceed as follows. First, to avoid dealing with the normalization of the scaling variable 
appearing in eqs. (EI), ([5]), and (EJ) we shall consider the unnormalized scaling function 
9(x, p) defined as 

9(tL 1 ^,p)=9(a tL 1 ^,p), (87) 

where in eq. (J3]) we have considered only the leading term in the expansion of u t ; here, 
according to eq. (EJ), ao is an amplitude not yet specified. 

We have already found in section H] that for p — > the amplitude of the critical 
Casimir force at T c exhibits a linear dependence in p. This leads us to introduce the 
following corresponding generalizations: 

9(T,p)=9(r,0) + pE(r)+O(p 2 ), (88) 

9(f,p) = 9(f,0)+pE(f) + O(p 2 ), (89) 

omitting possible higher-order terms in the preceding expansion. For each given pair 
(/3,L), inserting eq. ( 1881) or eq. ( 1891) into eq. ( 1691) leads to a linear dependence on p of 
the quantity I(j3,L,L\\). Accordingly, a simple linear interpolation scheme for the MC 
data from systems with different aspect ratios p yields, for every given pair (/3, L), 

x -,x 9(t(L- 1/2 + c) 1 / u ,p) , x 

1 09, L, p) = B(t) + K \ L _[ /2 + c)3 , (9°) 

with c from eq. (184")) . There are two reasons for interpolating the data at aspect ratio 
p instead of directly using the data at p as it was done in section HJ First, in this way 
we have suppressed the 1/L correction which arises from eq. (1701) . Secondly, as it will 
become clear in the following, by considering data for various lattice sizes with the same 
aspect ratio p, we are left with a dependence on L only in the first argument of the 
scaling function 9(t, p), a fact which enables us to eliminate the correction to scaling in 
an easy way. 

Next, we define a function g by taking the difference, at the same aspect ratio p, of 
the two expressions for / corresponding to lattice sizes L and aL, respectively: 

g{P, L, p) = (L- 1/2 + c) 3 \l 03, L, p) - I ((3, aL, p)} . (91) 
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As in eq. (1751) for figure EH we choose a = 2. Inserting eq. (19"U|) into eq. (19~T|) leads to 
g((3, L,p) = 9 (t(L - 1/2 + cfl\ p) - A(a, L)9 (B(a, L)t(L - 1/2 + c) 1 ^, p) , 
a, r^ f L-l/2 + c\ 3 D . r . faL - 1/2 + c\ l/v 

This eliminates the background term B(t). With a slight abuse of notation on the left 
hand side, we can rewrite the previous equation as 

g(f,p;L,a)=9(r,p)-A(a,L)9(B(a,L)f,p), f = t(L - 1/2 + c) l/u . (93) 

If we had simply used I instead of I, the two terms in eq. ( l93l would refer to the scaling 
function 9 calculated at different aspect ratios, leaving us with a complicated expression. 
With our choice a = 2, we have A(a, L) = l/8+0(l/L) and B(a, L) = 2 l ' u + 0(l/L) ~ 
3 + 0(1/ L). Keeping in mind that for large r the Casimir force decays exponentially, 
we see that in eq. (1931) the second term represents a correction to the first term, which 
vanishes in the limit a — > oo. In order to eliminate this correction we introduce 



g n (r, p; L, a) = J^ [A(a, L)f g([B(a, L)] k r, p; L, a), 

g (f, p; L, a) = g(f, p; L, a). (94) 

Inserting eq. (j93l) into eq. (1941 we obtain 

g n (f, p; L,a) = 6 (f , p) - [A(a, L)] n+1 6 ([B(a, L)} n+1 f , p) . (95) 

Thus we have g n (f,p;L,a) n ^°° t Q(f,p). Moreover, the error due to truncating the 
sum in eq. (l94"j) is proportional to [A(a = 2,L)] n ~ 2 _3 ( n+1 ). Accordingly, the sum 
in eq. ( 194]) converges quickly to the scaling function 9. Starting from a Monte Carlo 
estimate of g at a certain value of f , in order to be able to calculate the sum appearing 
in eq. (|94l) the value of the function g at f' — B k f is required, which might be not 
directly available from the MC data. However, this value can be estimated by using a 
simple interpolation spline for those values of the function g which are available. With 
the present precision of our data, at n = 2 the error associated with the truncation 
and given by eq. (|95l) is smaller than the statistical error bars. Therefore we use the 
approximation 9(f,p) ~ g2(f,p\L,a). We note that, if the procedure is correct, the 
reconstructed function 9 should not depend explicitly on a and L, but only on the scaling 
variable f and the aspect ratio p. Finally, we implement the appropriate normalization 
by using eq. ( IH71) with the normalization constant a^ given in eqs. (jSJ) and (15"I1) . 

In figures [5], [7J and we show the universal scaling function 9(t, p) for aspect 
ratios p = 1/6, 1/8, 1/10, and 1/12, respectively. In order to treat correctly the 
statistical covariance between the various quantities in eq. (l94"j) . the error bars have 
been calculated according to the jackknife procedure (see, e.g., Ref. [90]). We observe 
a good data collapse for L > 12, which supports qualitatively our procedure and the 
scaling ansatz in eq. ( 169]) . The data for L = 8 appear to be slightly off the curves 
obtained for larger lattices. This is not surprising because, as we already noted in 
section HI data for such a small lattice size suffer from higher-order scaling corrections. 
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Figure 5. The universal scaling function 8(t, p) of the critical Casimir force (eq. (|55jl ) 
for r = t{L/t^) l / v and for an aspect ratio p = 1/6. The error bars are smaller than 
the symbol sizes. 
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Figure 6. Same as figure [5] for the aspect ratio p = 1/ 
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Figure 7. Same as figure [5] for the aspect ratio p = 1/10. 
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Figure 8. Same as figure [S] for the aspect ratio p — 1/12. 
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Figure 9. The universal scaling function 9(x, p) of the critical Casimir force (eq. ([26])) 
as a function of the scaling variable r = i(X/£ c J~) 1 / 1 ' and extrapolated to the aspect 
ratio p — 0. We compare our results with the mean value (mv) of the critical Casimir 
forces for laterally homogeneous ++ and -\ — b.c, as obtained from the MC data 
in Ref. [60] for the three approximants I, II, and IV presented therein. There is 
satisfactory agreement with approximant IV, which was claimed to be the best one 
among those studied in Ref. [60] . 



In figure [9] we present the scaling function for the critical Casimir force extrapolated 
to p — > 0. As we mentioned at the end of section [3] in this limit the Casimir force is 
expected to be the mean value of the force for the laterally homogeneous ++ and H — 
configurations. The scaling functions for these b.c. have been computed by Monte Carlo 
simulations in Ref. [60], where 3 curves, denoted as I, II, and IV, have been presented. 
From those curves we have formed their mean values by using an interpolation spline 
for the three approximants presented in Ref. [60]. In figure [9] we show a comparison 
with our results. There is good agreement between our curve and the one obtained from 
the approximant IV of Ref. [60]. Incidentally, it is reported in Ref. [60] that among the 
three it is this approximant which describes the finite-size scaling of the MC data best. 

In figure [10] we show a comparison of the scaling function of the critical Casimir 
force for the aspect ratios p = 1/6, 1/8, 1/10, 1/12, 0, as obtained for L = 12, as well 
as the mean value of the force for the laterally homogeneous ++ and H — b.c. following 
from approximant IV of Ref. [60]. Within the available range of r values, the function 
9(t, p) increases for decreasing p. 
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Figure 10. Comparison of the universal scaling function 6(x, p) of the critical Casimir 
force (eq. ([26])) as a function of the scaling variable r = t(L/£o) 1 ' 1 ' for the aspect ratios 
p = 1/6, 1/8, 1/10, 1/12, and p — >■ for fixed L = 12. We also show a comparison 
of our results with the mean value (mv) of the critical Casimir forces for laterally 
homogeneous ++ and H — b.c, as obtained from the MC data in Ref. [60] for the 
approximant IV presented therein. 



5.2. Chemical- steps contribution 

As discussed at the end of section [3j in the limit p — > the critical Casimir force in the 
presence of a pair of individual chemical steps on one of the confining walls (figure [1]) 
reduces to the mean value of the forces for the laterally homogeneous ++ and H — 
configurations. This mean value can be interpreted as the force for an ersatz system in 
which the system shown in figure [1] is decomposed into two halves disconnected along 
the chemical steps which become infinitely separated in the limit p — > 0. Accordingly 
one expects that the effect of the presence of the pair of individual chemical steps on the 
critical Casimir force enters into its dependence on the aspect ratio p. This consideration 
can be formalized by generalizing the discussion presented at the end of section 12.31 In 
the following we first discuss how to define a line contribution to the singular part of the 
free energy density for general b.c, specializing later the argument to the present case. 
To this end we consider two slabs L x L\\ x L\\, which undergo a second-order phase 
transition, with the same bulk and surface universality classes. The confining surfaces 
of the first slab (1) exhibit laterally homogeneous b.c, such that the system displays 
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two surfaces but no edges. This can be realized by imposing periodic b.c. in both lateral 
directions and homogeneous b.c. on the two confining surfaces. The second slab (2) has 
one or more edges of extension L\\ or has a pair of individual inhomogeneities of linear 
extension L\\. The first case can be realized by imposing periodic b.c. only in one of the 
two lateral directions corresponding to the linear size Ly, leaving open b.c. in the other 
directions: such a slab has two surfaces of area Lii x L» and four edges of extension L». 
The second case can be realized by imposing a laterally inhomogeneous b.c. with linear 
spatial extent L», as in the system shown in figured! For these systems (1) and (2) 
away from criticality, i.e., for L>(, and in absence of external field H the free energy 
density T per volume LL'i and per fc^T decomposes as (compare eq. (jHJ)) 

F {1 \t,H = 0,L,L||) = / bulk (t) + i/surfft) + 0(e- L/ */L), 

J**\t,H = 0,L,L||) = / bulk (t) + y suTi (t) + j^hUt) + 0(e- L 't/L), (96) 

where p = L/L\\ and /bulk, /surf, and /i ine being independent of L and p. In the limit 
p — >■ the indicated correction terms in J 7 ^ and J-^ 2 ' reduce to the transversal finite 
size contributions 0(e~ L '^/L) which give rise to the critical Casimir forces in the lateral 
thermodynamic limit L\\ — > oo. For large but finite values of L\\ these correction terms 
acquire a dependence on p which depends on the b.c. of the system. 

Generalizing the arguments given in Ref. [76] and the corresponding discussion in 
section 12.31 we define the following quantity: 

/ line (t,L)^L 2 f|- [j**\t, H = U,L,L{)- 7**^1, H = Q,L,L{j\) , (97) 

\dpL,t l J\p=0 

so that, we have (compare with eq. ( I2TT) ) 

/ linc (t, L) = / linc (t) + 0(Le- L ^), L^oo. (98) 

Equation (|97|) formally defines a line free energy density /ij ne also in the critical regime. If 
the slab (1) is realized as in the example given above, the free energy of this slab contains 
only bulk, surface, and finite size corrections which have at most a quadratic dependence 
on the aspect ratio p; in such a case the definition in eq. (j9"7j) yields dJ 7 ^ 1 ^ /dp\ p=0 = 
and /ii n e can be identified as the only free energy contribution of the two systems which 
varies oc p for small aspect ratios (see also the discussion below). 

In the present case, we compare the free energy density for the geometry of a pair of 
individual chemical steps with the mean value of the free energy densities of the systems 
with laterally homogeneous b.c. ++ and H — . For these latter b.c. in the critical region 
the free energy density decomposes into (compare eq. (J2]) with / = / ++ or / = /+_) 

T++(t,H = 0,L,p) = fiZl(t) + y%$(t) + -^/ + +M,P), 

T+-(t,H = 0,L,p) = fl™l(t) + i/2(t) + -^/+-M,p), 

r = UtL 1 ^, (99) 
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where /++(r, 0, p) and /+_(t, 0, p) are scaling functions which describe the singular part 
of the free energy density for the laterally homogeneous ++ and H — b.c, respectively, 
in the absence of a bulk fielcp. As in eq. ( }9T|) . the free energy density / ste ps for the pair 
of chemical steps is defined by 

° F(t,H = 0,L,p) 



J steps \ti L) 



dp 



L.t 



-±J*+(t,H = 0,L,p)-±J*~{t,H = 0,L,p) 



(100) 



p=0 



with J-(t, L, p) as the free energy density for the system shown in figure [U Its singular 

' steps 



Part /steps is S iven °y ( see eq. &) 



fiU*' 1 )- L { dp 



/(r,0,p)--/ 4 



(r,0,p)-i/+_(r,0,p) 



(101) 



p=0 



where /(r, 0, p) is the scaling function of the singular part of the free energy density 
for the system in the pair of individual chemical steps geometry and in absence of 
bulk field (see eq. (J2j)). In eqs. ( 11001) and ( 11011) we indicate with ++ and H — the 
quantities relative to the systems with laterally homogeneous ++ and H — b.c, while 
conforming to the notation of the previous sections the corresponding quantities for the 
system shown in figure CD are indicated without further specifications. Equation f II 1 ft 
renders the following relation, in lowest order in p, between the three functions /(r, 0, p), 
f++(r,0,p), and / + _(r,0,p): 

/(r,0,p) = i/ ++ (r,0,p) + i/+_(r,0,p) + P Lf% s (t,L) + 0(p 2 ), p -► 0,T,L fixed. 

(102) 

Note that Lff eps (t, L) L »g, Lf^l ps (t) ~ L|i| 2 - Q < ~ |t|", due to a, = a + 2i/ [68J, so that 
the r.h.s. of eq. f ll02[) is indeed a function of p and r. Together with the definition of 
the critical Casimir force in eq. (1231) one has 
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dL 
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where 
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c,- 
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<9L 




t,L || 


9((l/L 2 )/ ++ (r,0,p)- 


- £/&(*)) 


<9L 




a((l/L 2 )/ + _(r,0,p)- 


" */&(*)) 



<9L 



i,Ln 



i,Ln 



(103) 

(104) 
(105) 
(106) 



t) We note that in the absence of a bulk field, the surface free energy densities are the same for both 
b.c. 
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are the Casimir forces in the pair of individual chemical steps, ++, and H — geometry, 
respectively. The last term in eq. (I103P is the contribution insteps to the critical Casimir 
force due to the pair of individual chemical steps: 



d[(p/L)fiL(t,L) 



C, steps 



dL 



t,L\ 



+ 0( " 2) --u 



dfstepaV'iL) 



OL 



+ 0(p 2 ). (107) 



t,Lii 



According to the Monte Carlo results in Ref . [60] , the dependence on p of the critical 
Casimir force for the laterally homogeneous ++ b.c. is negligible for p < 1/6, while for 
the H — b.c. it is quadratic and becomes relevant in the low-temperature phase. (Note 
that this quadratic term has a zero at r = 0, see eq. f l30|) .) This implies that eq. f llOip 
reduces to 

l<9/(r,0,p) 



Kb) 

J steps 



(t,L) 



L 



dp 



(108) 



p=0 



Inserting this result into eq. (11071) we obtain 
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On the other hand from Fq 
one finds from eq. (|26|) 



so that 
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T ddf(T,0,p) 
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v dr dp 
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The comparison with eq. (1109J) yields 

p de{r,p) 



C,steps 



L 3 dp 



+ 0(p 2 



p=0 



With the expansion in eq. (188|) we obtain 
F c , steps = ^E(r)+0(p 2 ), 



r 



t{L/itf lv - 



:in^ 



(112) 



(113) 



Since, as stated above, the dependence of (Fc,++ + i 7 c ]+ _)/2 on p is quadratic, eq. f ll 13j) 
implies that in the limit p — > the contribution to Fc, which is linear in p, is solely due 
to the presence of the pair of individual chemical steps on one of the confining surfaces 
and thus serves as its fingerprint on the critical Casimir force. 

We can extract this contribution from the MC data. From eqs. ( 169|) . ( 170|) . and (IHHj) 
we have 



dI(J3,L,L\\ 
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l/v^ 



(114) 
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Figure 11. The scaling function E(t = t{L/(^) 1 / v ) which describes the contribution 
of the pair of individual chemical steps to the critical Casimir force via the dependence 
of the latter on the aspect ratio p = L/Ln (see eqs. (|88|) . (|103[) . and (|113[1 ). The 
statistical error bars for L — 8 and L = 12 are smaller than the symbol size. 



We note that the coupling parameter approach outlined in section [3] results in the free 
energy difference between two systems with the same chemical steps. Therefore in 
I(J3,L,L\\) non-singular background terms in the surface and line contributions to the 
free energy drop out. Since for every pair /3, L we have simulated systems of various 
aspect ratios, the derivative in eq. (j!14p can be inferred from a simple fit linear in p, 
which is the same fit as the one which has been used in order to process the quantity 
I(f3,L,L\\) in eq. ( 190]) . From this the function E{r) follows according to eq. (j!14p . 

In figure [Tl] we show the function E(t) as extracted from our MC data. It is 
negative within the whole range of r values explored, which is consistent with the fact 
that the critical Casimir force is a decreasing function upon increasing p. The data for 
the various lattice sizes displayed in figure HH collapse onto a single curve, with slight 
deviations for the data from the smallest lattice size. This is consistent with the findings 
of section H] and section I5.1[ according to which the data for L = 8 are influenced by 
subleading scaling corrections. Moreover the data collapse confirms that the quantity 
we have extracted as the contribution of the pair of individual chemical steps indeed 
represents the singular part only, described by the scaling function E(t) in eq. (I113p . 
Since the critical Casimir force decays to if \t\ — > oo, which holds for an arbitrary 
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aspect ratio p, it follows from eq. (18"8"]) that E(t — > ±00) = 0. This is confirmed in 
figure [11] for r — > +00, while the limiting behaviour for r — > —00 lies outside the 
available MC data. Therefore we expect that E(r) reaches a minimum at r = tq < —8. 

According to the reconstruction scheme (see eq. ([91)1 ). in order to extract the scaling 
function 9(r, p) at some value r ^0 the system has to be simulated at [B(a, L)] n r, with 
n given by the truncation of the sum in eq. f )94[) . In particular, for the present data one 
has n = 2 and B(a,L) ~ 3 (see section |5\T|) ; thus in order to calculate 9(r,p) within 
the range —8 < r < 0, one needs data for —72 < r < 0, where the data in the interval 
—72 < r < — 8 are used only to reconstruct the function 9(r, p). Although at sufficiently 
low temperatures the function 9(r, p) is suppressed, one still needs simulations for a large 
interval in the low-temperature phase. Since the computational cost of the simulations in 
the low-temperature phase increases with decreasing temperature, these circumstances 
limit the availability of data for T < T c . 

The film geometry studied here is relevant for the critical Casimir force in the 
presence of a chemically structured substrate. The simplest realization of such a 
substrate consists of a substrate which is finite in one lateral direction (x) and has 
a macroscopic extent in the other direction (y) so that it is de facto translationally 
invariant in this latter direction, exhibiting a single chemical step in x direction (see 
figured]). For such a substrate the film geometry can be approximately realized by either 
considering a wetting film of a binary liquid mixture p2][T3] or by a colloidal particle in 
front of such substrate [17J. In the first case the wetting film thicknesses forming next to 
a + or — surface adjust to the different corresponding substrate potentials, so that the 
resulting critical Casimir force is described by the film geometry considered here only 
if the concomitant non-uniformity of the film thickness is small. This is conceivable 
because the wetting film thickness is mainly determined by the total density whereas 
the critical Casimir force is linked to the concentration fluctuations. In the second case 
the film geometry is approximately recovered if the radius of the colloidal particle is 
much larger than its distance from the substrate. In both cases the critical Casimir 
force is influenced by the lateral b.c. which in a first approximation could be treated as 
open ones. As we already mentioned in section [1] for such a geometry the aspect-ratio 
dependence of the critical Casimir force is due to the presence of the chemical step 
and of the lateral edges, the contributions of which in general cannot be disentagled. 
On the other hand, the linear aspect-ratio dependence of the critical Casimir force 
calculated in this section is due to the presence of two individual chemical steps, while 
the corresponding contribution due to lateral edges could be determined in a similar way 
by considering a film geometry in the presence of homogeneous surfaces and laterally 
open b.c, or b.c. adapted to the actual experimental conditions realized there. Then 
the expected critical Casimir force Fq for the film geometry in the presence of a single 
chemical step and laterally open b.c. (or b.c. adapted to the actual experimental 
conditions) is given by 
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c L 3 



9(t, 0) + p£ edges (T) + ^E(t) + p6E(rj\ + 0(p 2 ), 



where E e( x ges (T) is the scaling function associated with the aspect-ratio dependence due 
to the lateral edges, to be determined as described above, 9(t, p) and E{r) are the scaling 
funtions calculated here, and 5E(t) is a scaling function which accounts for the expected 
non-additivity of the aspect-ratio dependence of the critical Casimir force with respect 
to the individual line contributions. One might expect that SE(t) is small compared 
with the other three contributions. 

6. Order parameter profiles 

The free energy and the critical Casimir forces provide integral informations about finite 
sized systems. Order parameter distributions deliver valuable additional and spatially 
resolved informations, which provide a deeper understanding and predictions which can 
be probed experimentally, e.g., by X-ray scattering under grazing incidence (see, e.g., 
Ref. [91]). 

We have computed the order parameter profiles for the system shown in figure [1] 
at the critical temperature, for lattice sizes L = 16, 24, 32, 40, 48 and aspect ratios 
p = 1/6, 1/8, 1/10, 1/12. We have sampled the order parameter in a region close to the 
"central" chemical step located at z = x = (see section E]), so that the presence of the 
"lateral" chemical step induced by the lateral periodicity is not relevant for the results 
here. In the following by referring to a chemical step we mean the central one. Certain 



details of the simulations are reported in Appendix B The leading scaling behaviour of 



the order parameter profiles $(x,z,L,L\\) at criticality can be obtained from eq. ([36 

by setting t = 0: 



<!>(x,z,L,L ll ) = B\- I \ «^--p), (116) 

with (see eq. ( 13~7|) ) 

cj) c (x, z, p) = <f>± (x, z, 0, p) , (117) 

where we have introduced the scaling function C (x, z, p) in order to simplify the 
notation. In addition to the leading scaling behaviour as of eq. (I116p . we observe 
scaling corrections oc 1/L. However, within the precision of our data, no aspect ratio 
dependence has been found in the spatial region near the chemical step which we have 
considered here. In order to extract the thermodynamic limit, we fit our MC data for 
$(x, z, L, L||) to the following expression: 

HxM ^ B (±y'"^,±) + i g ^)y ( 118 ) 

where, with a slight abuse of notation, we have dropped the dependence on the aspect 
ratio p, and (f) c and g are unconstrained functions. This is achieved by performing a 
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fit linear in 1/L of $(x, z, L)L^I V for every pair (x/L, z/L), using (3/u = (1 + r/)/2 = 
0.51819(8) [71]. Error bars have been determined from the jackknife procedure (see, e.g., 
Ref. [90]) in order to take into account the statistical covariance of the data sampled from 
the same MC run for various pairs (x/L, z/L). The non- universal amplitude combination 
appearing in eq. (11 18ft is given by (see eqs. (1M1) and (135]) ) 

B(£+f /v = 0.933(6) (119) 

and has been computed by analyzing the 25 th -order high-temperature expansion 



reported in Ref. [7T] (see Appendix A for details). In sections H] and [5] we observed 
that data for the lattice size L = 8 are affected by next-to-leading scaling corrections, 
while for L > 10 no subleading scaling corrections are observed. In the present data for 
the order parameter profiles the minimum lattice size is L = 16, so that corrections to 
eq. (I118P should be negligible. In fact, our fits always render a good \ 2 1 DOF, except 
for a narrow spatial region close to the chemical step. This indicates that there effects 
due to subleading corrections are probably strong. In fact, if scaling corrections are 
(at least partially) due to analytic corrections of the scaling variables as discussed in 
section [21 the function g(x/L, z/L) should be expressable in terms of derivatives of the 
function <f) c (x/L,z/L) (compare with eq. (133]) ). Since close to the chemical step the 
function <p c (x/L,z/L) varies steeply, the ensuing leading correction to scaling oc L^ 1 
can potentially become large enough as to invalidate the ansatz given in eq. (11191) . 

In figure [12] we report the function <j) c (x/L, z/L), for various values of z/L in the 
region close to the chemical step located at x — 0. For our simulations we have chosen 
the coordinate system such that the confining surface with the spins fixed to form a 
chemical step corresponds to z = 0, while the homogeneous surface with the spins fixed 
to +1 is located at z = L — 1. In the lateral direction the origin x = corresponds, in 
the lower surface, to spins fixed to +1, whose left neighbours are spins fixed to — l lffl . 

By inspecting the profiles we see that close to the chemical step at \x/L\ <C 1 and 
z/L = 1/8, the function varies strongly, passing from negative values at x/L < to 
positive values at x/L > 0. As expected, upon increasing z/L the gradient at x/L = 
decreases, as the frustration due to the chemical step heals. At z/L = 7/8, i.e., close to 
the laterally homogeneous boundary, the scaling function is almost flat, signalling that 
there the effect due to the presence of the chemical step on the distant wall is small. We 
also note that for \x/L\ > 0.5, there is de facto no dependence of the order parameter 
profiles on the lateral coordinate x: at such distances from the chemical step, the order 
parameter is no longer influenced by the chemical step and we recover the profiles (as a 
function of the coordinate z normal to the surface) for homogeneous boundaries: — \- for 
x/L < —0.5 and ++ for x/L > 0.5. Consistently, for x/L > 0.5 the scaling function is 
antisymmetric around z/L = 1/2 and for x/L > 0.5 it is symmetric around z/L = 1/2. 
These findings support the line of reasoning at the end of section [3] concerning the 

ff In the finite size scaling limit, i.e., in the limit L — > oo at fixed z/L, x/L, the precise choice of the 
coordinate system (e.g., whether the lower boundary is at z = or z = 1) does not affect <fi c (x/L, z/L); 
it only influences the correction-to-scaling function g (x/L, z/L). 
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Figure 12. The universal order parameter scaling function <f> c (x/L, z/L) at criticality 
as denned by eq. (|116[) . For x < (x > 0) there are — (+) b.c. on the lower confining 
surface at z — whereas there are + b.c. on the upper surface at z = L — 1. The error 
bars are smaller than the symbol size. The data points are interpolated with a cubic 
spline. 



various contributions to the critical Casimir force, even at T c . 

Close to the walls and in the regions where the profiles resemble those for laterally 
homogeneous b.c, the scaling of the order parameter is predicted to exhibit the distant- 
wall corrections described by eqs. ( 1461) . ( H71) . and ( )48l) . We first consider the case of 
++ b.c. To this end we analyze the data close to the bottom wall and for x/L > 1 in 
order to avoid a potential bias due to a residual dependence on x/L. We fit them to 
the expression 



${x,z,L) = CL^ jh 



-PI" 



1 + A 



1 + j9(z/L) 



(120) 



with C, A, and g(z/L) as free parameters. The correction term ~ {z/L) 3 is expected to 
be valid only if z/L is sufficiently small. Thus we have performed the fits by considering 
only those data with (z/L) < (z/L) m ^. The corresponding results are reported in 
table SJ 

By analyzing conservatively the dependence of the results of the fit on (z/L) max , 
we obtain the following estimates: 



C = 0.7879(5), 



121^ 



Critical Casimir forces and adsorption profiles 



38 



Table 4. Fits of the order parameter profile in the region with ++ b.c., i.e., x/L > 1 
and (z/L) < (z/L) max to eq. (|120|) . The indicated error bars are the sum of the 
statistical uncertainty (first number, obtained with the jackknife technique [90| ) and 
of the variation of (3/v = (1 + i])/2 = 0.51819(8) [71] due to its error bar (second 
number). 



(•Z/ijmax 


C A 


g{z/L) 




X 2 /DOF 


3/8 


0.7876(2 + 1) 1.442(7 + 2) 


5(1/8) = 
.9(2/8) = 
5(3/8) = 


-1.825(3 + 2) 

-0.798(4+1) 
-0.338(8 + 2) 


1.24 



4/8 



0.7878(2 + 2) 1.408(3 + 1) 



5(V8) 

5(2/8) 
5(3/8) 
5(4/8) 



-1.831(4+2) 

-0.795(4+1) 

-0.309(5+1) 

0.060(9 + 2) 



1.26 



5/8 



0.7881(2 + 2) 1.393(2 + 1) 



5(1/8) 
5(2/8) 
5(3/8) 
5(4/8) 
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-1.837(4- 
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0.089(5 - 
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1) 
1) 
2) 
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1.24 



6/8 



0.7878(2 + 2) 1.4015(7 + 3) 



5(V8) 

5(2/8) 
5(3/8) 
5(4/8) 
5(5/8) 
5(6/8) 



-1.831(4- 

-0.792(4- 

-0.302(4- 

0.076(4 - 

0.423(4- 

1.019(6- 



1.26 



7/8 



0.7810(1) 1.5146(3) 



5(1/8) 
5(2/8) 
5(3/8) 
5(4/8) 
5(5/8) 
5(6/8) 
5(7/8) 



-1.662(4) 
-0.654(4) 
-0.247(4) 
-0.031(4) 

0.051(3) 
0.267(2) 
2.268(2) 



56.7 



A= 1.40(1), 
(7(1/8) = -1.834(8), 
g(2/8) = -0.794(6), 
3(3/8) = -0.305(9), 
^(4/8) = 0.07(2), 
0(5/8) = 0.44(2). 



(122) 
(123) 
(124) 
(125) 
(126) 
(127) 



By comparing eq. ffT20|) with eqs. (jSJ, Q3BJ, and (gSJ one has C = B(^f/ U c + . By 
inserting the estimates given in eqs. ( II 19ft and ( 1121J) we obtain 

c+ = 0.844(6). (128) 

We can compare this result with previous theoretical estimates c + = 0.87(7) |92j, 
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c + = 0.94(5) |H2|, and c + = 0.857 [03] • From a variety of experimental data 
corresponding to different liquids one infers c + = 1.60(42), c + = 0.77(19), c + = 1.14(29), 
c+ = 0.91(26), c+ = 1.05(9), c+ = 1.02(10), c+ = 1.25(9), c+ = 0.84(15) [HS]. Our 
result is one order of magnitude more precise than previous theoretical determinations 
and in very good agreement with the results in Refs. [921J93] , while the discrepancy with 
Ref. [82J is of only two error bars. The comparison with the experimental data extracted 
in Ref. [82] is less satisfactory: half of the experimental data are in agreement with our 
result, while the apparently more precise determinations are not compatible with our 
result. 

The comparison of eq. (11201) with eqs. (J47J), (HE]), and (BED yields 
A = — [d — 1)A ++ C + , where A ++ is the critical Casimir amplitude for ++ b.c, and 
C + is a universal coefficient which depends only on the surface universality class of the 
close surface, i.e., here the extraordinary or normal surface universality class. With 
{d- 1)A++ = -0.820(15) [64] we obtain 

C + = 1.71(4). (129) 

This coefficient C + has been determined for the extraordinary surface universality 
class in Ref. [84] , using a combination of the e-expansion technique and interpolation 
with the exact d = 2 result leading to C + ~ 1.42 — 1.96. Our determination is more 
precise, and in full agreement with this estimate. 

We have repeated the above procedure for the data of the profile for x/L < —1, 
corresponding to H — b.c. In this case the fits of eq. (I120p result in a large y 2 / DOF. 
Due to (d — l)A+_ = 5.613(20) [64], in this case the distant wall gives rise to a correction 
which is about seven times stronger that in the ++ case. Using the estimate in eq. ( 11291) 
and the Casimir amplitudes (d- 1)A++ = -0.820(15) and (d- l)A+_ = 5.613(20) [M], 
we can infer that eq. (T4H1) is presumably valid for distances z such that 

C + (d-l)A +± (£) 3 «l, (130) 

y<0.9, for + + b.c, (131) 

Lj 

y<0.5, for + - b.c. (132) 

Equation H132J) explains why fits to eq. f 1 1 2 j) fail for -\ — b.c. for the presently available 
MC data. In fact, only the data for z/L = 1/8 and 2/8 can satisfy eq. (1132|) . but they 
are insufficient for performing a fit to eq. (I120p . On the other hand, the bound given by 
eq. (I13ip is consistent with the results in table HI where fits to eq. f 1 1 2 j) turn out to be 
reliable up to (z/L) max = 6/8 = 0.75. 

7. Summary and conclusions 

We have studied the critical Casimir force and the order parameter profiles for a three- 
dimensional slab of thickness L belonging to the Ising universality class. We have 
employed periodic boundary conditions in the two lateral directions of extents L\\ and 
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fixed boundary conditions on the two confining surfaces. The Ising spins on the upper 
surface are fixed to +1. The lower surface is divided into two halves, one with spins 
fixed to —1 and the other with spins fixed to +1 (see figured]). We have investigated 
this system by combining Monte Carlo (MC) simulations and finite-size scaling analysis. 
We have obtained the following main results: 

• In the finite-size scaling limit t = (T — T c )/T c — > 0, L — > oo, with £/L and L/L» 
fixed, the singular part of the critical Casimir force Fc per area L 2 and per k^T 
takes on the universal scaling form Fc(t, L, L\\) = L~ 3 9(t = t(L/^Y^ u , p = L/L\\) 
(eq. f j26|) ) with £(£ -> O 1 * 1 ) = c^li] -1 ' as the true bulk correlation length and p as 
the aspect ratio of the slab. The force Fc is determined by integrating the thermal 
Monte Carlo average of a suitable crossover Hamiltonian (see equations (!55j) . (|68|) . 
and (1691)). 

• In order to be able to extract the universal scaling function 9(t, p) from the MC data 
for a finite size system, particular care has been taken to minimize the influence of 
corrections to scaling. To this end we have studied a suitable Blume-Capel model 
which belongs to the Ising universality class and which suppresses the leading 
correction to scaling oc L -0 - 832 (eq. (j4~9l ). From the available high-temperature 
series of this model we have determined the relevant non-universal amplitudes 



( Appendix A ). This allows one to properly normalize the scaling functions and 
scaling variables. A detailed finite-size scaling analysis leads to a relation between 
the actual MC observable called / (eq. (155]) ) and the desired scaling function 9(r, p) 
with suitably adjusted scaling variables in order to minimize corrections to scaling 
(eq. (!69l). This procedure yields numerically accurate data for the scaling function 
0(r,p). 

As expected on general grounds (figure [2]), in the limit of vanishing aspect ratio 
p = L/L\\ —¥ the critical Casimir force for the system shown in figure [I] reduces 
to the mean value of the critical Casimir forces for laterally homogeneous ++ and 
H — boundary conditions so that 9(r,p — > 0) = ~[0 ++ (r) + + _(r)] (see figure [9] 
and p — > and mv(IV) in figure [TUD. Since Fc t +- is more repulsive than Fc,++ is 
attractive, 9(t, p — > 0) is positive (figures 191 and fTDl so that for p — >■ the critical 
Casimir force for the system in figured] is repulsive. 

The presence of the pair of individual chemical steps documents itself in the 
dependence of the scaling function 9(r,p) on the aspect ratio p. For laterally 
homogeneous b.c. in the limit p — > the critical Casimir force generally approaches 
its limiting value oc p 2 . The amplitude of this quadratic term vanishes at T c 
(see eq. (130])). The presence of the pair of individual chemical steps generates 
an additional line contribution oc L\\ to the free energy (eq. (196]) ) which causes 
a leading, linear dependence on p (figure E]). Based on a suitable extrapolation 
scheme, at T c this linear variation of 9{t = 0,p) = 0(p) is shown in figure H] 
Since for p — > the linear dependence on p dominates, the decrease of 0(p) upon 
increasing p implies that the presence of the chemical steps weakens the repulsive 
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critical Casimir force relative to the mean value of the forces for the corresponding 
laterally homogeneous b.c. 

• For T 7^ T c we have determined the scaling function 0(r,p) for p = 1/6, 1/8, 1/10, 
1/12 (figure El El d and El respectively). It turns out that 9(r, p) < 9(r,p = 0) 
for the available values of the scaling variables r with the deviation being most 
pronounced for r < (figure [TUT) . This difference is captured by the scaling function 
E(t) defined via Q[r,p ->• 0) = 6{r,p = 0) + pE{r) + 0{p 2 ) with E(t) < and 
E(t — ¥ ±oo) = (see figure [TTJ and the discussion at the end of section IST2"]) . This 
linear contribution ~ E(t) is solely due to the pair of individual chemical steps. 
Also for this laterally inhomogeneous system the amplitude of the quadratic term 
vanishes at T c (see eq. (1501)). 

• At the critical temperature we have determined the order parameter profiles for 
the system shown in figure [TJ (see figure [12]) . Since in the spatial region around 
the chemical step which we have considered here the dependence of the order 
parameter profiles on the aspect ratio is weak, these profiles can be described by 
the scaling function (f> c (x/L,z/L, p = 0) (see eq. (I116P ). Sufficiently away from the 
chemical step the order parameter profiles reduce to the ones which correspond to 
the laterally homogeneous b.c. ++ or H — . For the latter ones we have determined 
the universal amplitude c + (equation (11281) ) characterizing the leading behaviour 
of the order parameter near a wall (equation (j4"6l) ) for which experimental data 
are available [82J, as well as the amplitude of the leading distant wall correction 
(equations (1171) and (14"8~|) . and ( 11291) ). which is in full agreement with previous, but 
less accurate, estimates. 

Our results are relevant for the critical behaviour of confined systems belonging 
to the Ising universality class and in the presence of a chemically structured substrate. 
Indeed, such a system has been experimentally realized; the critical Casimir force has 
been probed by a spherical colloidal particle close to a structured substrate [TT|IT8"]. 
The film geometry studied here is realized approximately by such a system if the 
radius of the colloid is large compared to its distance from the substrate [53] . Another 
possibility would be to monitor the thickness of a wetting film of a classical binary 
liquid mixture near its critical end point of demixing in equilibrium with its vapour 
phase [T5] and in contact with a chemically structured substrate. However, in such 
a system in addition the wetting film thicknesses forming next to a + or — surface 
adjust to different adsorption preference. This adds a new interesting aspect to the 
problem. Some additional considerations regarding a possible experimental realization 
are reported at the end of section 15.21 

The present study points towards several interesting issues to be investigated in the 
future. A natural generalization consists of studying the critical Casimir force in the 
presence of many stripes with alternating adsorption preferences. In this case the critical 
Casimir force depends additionally on the width of the stripes (see Ref. [53] where the 
critical Casimir force between such a substrate and a colloid has been studied; genuine 
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three-dimensional simulation data for this system are still missing). Furthermore, one 
can consider an alternating adsorption preference also on the upper boundary, in which 
case even for a slab also a lateral Casimir force arises. Since such a system has been 
studied within mean field theory in Ref. [H], a comparison with three-dimensional 
simulation data would be of particular interest. The critical Casimir forces for more 
complicated geometries, such as the square-patterned substrate experimentally realized 
in Ref. [17], have not yet been studied theoretically. 

From a more theoretical point of view, the nature and the origin of the scaling 
corrections for non-periodic boundary conditions calls for further analysis. In particular, 
so far the ansatz of eq. (1321) appears to describe correctly those additional corrections 
to scaling proportional to 1/L which emerge in the presence of non-periodic boundary 
conditions. This has been also checked numerically in Refs. [6i~ | l80 l l8T] for the XY model 
with free surfaces and in Ref. [M] for the Ising universality class with fixed surfaces. A 
deeper theoretical understanding of this fact in terms of renormalization group theory 
would be highly welcome. 
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Appendix A. Analysis of the high-temperature series 

In this appendix we analyze the 25 th -order high-temperature expansion for the improved 
Blume-Capel model reported in Ref. [7J], with the aim of calculating non-universal 
amplitudes associated with various observables, which allows us to properly normalize 
scaling variables and universal scaling functions. We follow closely the notation in 
Refs. [67J[7I]. We consider the magnetization M per volume, the two-point correlation 
function G(x) and the corresponding moments m 2 j, and the four-point susceptibility Xi 
defined as 

M^£>(aO>, (A.l) 

X 

G(x) = (S(0)S(x)}, (A.2) 

m 2j = Y,\x\ 2j G{x), (A.3) 

X 

X4 = J2 (S(0)S(x 1 )S(x 2 )S(x 3 )) c , (A.4) 

2:1,22,2:3 
where the subscript c in the definition of %4 indicates the connected part of the thermal 
average of the product of spins. The susceptibility x an d the second-moment correlation 
length £ 2 nd follow from 

X = m , (A.5) 
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and 

&d = g- (A-6) 

Close to the critical temperature these quantities exhibit the following singular 
behaviours: 

X =C + r\ t>0, (A.7) 

6nd = C&J-", t > 0, (A.8) 

X4 = -C 4 + £" 74 , £>0 (A.9) 

M =B(-t)~ , £<0, (A.10) 

where t = (T — T c )/T c — ^ is the reduced temperature. We note that the exponential 
(or true) correlation length £, which governs the exponential decay of G(x), diverges for 
t -> 0+ as 

£ = #r", £>0, (A.ll) 

where the ratio £o7£j 2 nd = 1-000200(3) [H] is universal. 

We analyze the high-temperature expansion using quasi-diagonal first- and second- 
order integral approximants (IAl and IA2, respectively) [HI]. IAl are the solutions f(x) 
of the first-order differential equation 

Pi(x)f(x) + Po(x)f(x) + R(x) = 0, (A.12) 

where P\(x), Pq(x), and R(x) are polynomials of order mi, mo, and k, respectively. The 
overall normalization of the equation is fixed by setting -Pi(O) = 1 and the coefficients 
of P\(x), Po{x), and R(x) are fixed by the requirement that the Taylor expansion of 
f(x = (3 — > 0) matches the known high-temperature expansion. If x \ is the smallest 
real positive root of P\[x\ for x — > x c i the solution of eq. (1A.12I) behaves as 

f(x) ~ Ai(x)\x -x c i| 61 + Bi(x), x ->■ x c i, 
Pq{x c i) 



ei 



A'^cl) 



5l(Xel) = -^-, (A.13) 

where, near x = x cl , Ai(x) and -Bi(x) are analytic functions which are determined by 
P {x), P x (x), and R{x). 

IA2 are solutions of the second-order differential equation 

P2(x)f"(x) + Pi(x)f(x) + P (x)f(x) + R(x) = 0, (A.14) 

where P2(x), Pi(x), Po(x), and R(x) are polynomials of order m 2 , mi, m , and k, 
respectively, and ^(0) = 1. Again, the coefficients of the polynomials in eq. (IA.14J) are 
determined by the high-temperature expansion of f(x). The solution of eq. ( 1A.14I) has 
a structure which is similar to the one in eq. (1A.13I) : close the smallest real positive root 



Critical Casimir forces and adsorption profiles 44 

x c2 of P 2 (x) one has 

f(x) ~ A 2 (x)\x - x c2 \ e2 + B 2 (x), x ->■ x c2 , 

Pi(x c2 ) 

e 2 : 



P^(x c2 ) 



b ^ - -IS- (A - i5 » 

Inspection of eqs. ( 1A.13I) and ( IA.15I) tells that the critical exponents e\ and e 2 and 
the background terms Bi(x c i) and B 2 (x c2 ) can be obtained without solving explicitly 
eqs. (IA.12I) and (1A.14I) . because it is sufficient to determine the roots of P\(x) and P2(x), 
respectively. 

However, we are interested in the non-universal amplitudes -A$(x), which cannot 
be extracted directly from the coefficients of the polynomials P%{x) and R(x). In order 
to obtain them, we use the method outlined in Ref. [95J. Given the high-temperature 
series for a quantity Q((3) which close to criticality diverges as Q(f3) — A\(3 c /(3 — l\~ a , 
we analyze the series of Q((3) = Q(/3)(l — (3/(3 c ) a . The background term obtained 
from the analysis of Q(j3) corresponds to the desired amplitude A{j5) of the quantity 
Q((3). This approach requires the knowledge of the critical inverse temperature (3 C as 
well as of the critical exponent a which characterizes the critical behaviour of the given 
quantity Q(f3). As in Ref. [71], we consider quasi- diagonal approximants IA1 and IA2, 
i.e., integral approximants constructed from polynomials Pi(x), R(x) of almost equal 
degree, which are expected to lead to more reliable results [94]. If the series expansion 
of Q(J3) is known up to the order /3 n , we limit the orders of the polynomials Pi(x) and 
R(x) such that in the case of IAl we have 

Max{(n-2)/3-g,2} <m 0) m l) k < (n-2)/3 + q, (A.16) 

and in the case of IA2 

Max{(n-4)/4-g,2} < m ,m 1 ,m 2 ,k < (n-4)/4 + g. (A. 17) 

In eqs. ( 1A.16I) and ( 1A.17I) q measures the off-diagonality allowed, i.e., it limits the 
difference between the degrees of the polynomials Pi(x) and R(x). We have considered 
q = 3 for IAl and q = 2 for IA2. In order to improve the reliability of the result, we 
have always considered biased approximants, i.e., we set P\(x) = (1 — x/f3 c )Pi(x) in 
eq. (1A.12J) and P 2 (x) = (1 — x/(3 c )P 2 (x) in eq. (1A.14[) . so that the singularity occurs at 



x = (3 C . Furthermore, as in Refs. [7IJ[96] we discard the approximants which lead to 
spurious singularities in the region of the complex plane given by 

a; m in < Re z < x max , |Im z\ < y max , (A. 18) 

where z = (3/(3 c and the size of the rectangle is chosen as in Ref. [96J as x m ; n = 0.5, 
x max = 1.5, and y max = 0.5, where this choice was motivated by considering stability 
criteria. 

In Ref. [7T] the high-temperature series of xW) f° r f ne improved Blume-Capel model 
is reported up to the 25 th order. With the method described above, we have determined 
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the non- universal amplitude C + appearing in eq. (1A.7J) . With IAl approximants we 
obtain 

C + = 0.465994(3) ± 3 x 10~ 4 ± 6 x 10 -4 , for IAl, (A.19) 

where the first quoted error indicates the spread of approximants (standard deviation), 
the second stems from the variation of (3 C = 0.3856717(10) [71] within one error bar, and 
the third one originates from the error bar of 7 = 1.2373(2) [71]. With IA2 approximants 
we obtain 

C + = 0.4664(7) ± 4 x 10~ 4 ± 6 x 10" 4 , for IA2. (A.20) 

As a final estimate we take C + = 0.466(2). 

The analysis of the 21 th -order series of xa reported in Ref. [71] gives the following 
results: 

C+ = 0.3683(7) ± 3 x 10" 4 ± 10" 3 , for IAl, (A.21) 

and 

C+ = 0.365(4) ± 6 x 10" 4 ± 10~ 3 , for IA2, (A.22) 

where again the first error indicates the spread of approximants, the second one is due 
to the uncertainty of j3 c , and the third one follows from the spread of 74 = 27 + 3z/ = 
4.3650(6) [71]. As a final estimate we take C^ = 0.365(5). This allows us to extract 
the amplitude £o2nd °f the second-moment correlation length appearing in eq. ( 1A.8J) by 
using the universal amplitude ratio g\ [6] 

ql_ 



'/■' = ,JL ,„ (A. 23) 

1,2nd J 



which corresponds to the critical value of the zero-momentum four-point coupling 
constant. From Ref. [71] we quote g~£ = 23.56(2). Using this value and the ones 
for the amplitudes C + and Cf we finally obtain 

&nd = 0.415(2). (A.24) 

Using the universal ratio 6oV6o2nd = 1-000200(3) [71], from this one can calculate 
the non-universal amplitude of the true correlation length which, given the available 
precision, is identical to eq. ( 1A.24J) . We have also analyzed directly the series for 
£ 2 //3 = m2{(3)/(/3x(P))- The corresponding result ^ 2 nd = 0.419(7) appears to be less 
precise, but is in agreement with eq. ( 1A.24J) . 

The non-universal amplitude B in eq. (IA.10I) can obtained by using the universal 
amplitude-ratio R± [67] : 



r+R 2 

From Ref. [71] we quote R\ = 7.81(2). Using this value leads to 

5 = 1.47(2). (A.26) 

Finally, from eqs. ( fA~24l) and (fA~26l) together with /3/u = (1 + r])/2 = 0.51819(8) [71] 



we obtain the coefficient appearing in eq. (II 18ft : 

B{Z+Y' V = 0.933(6). (A.27) 
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Table Bl. The total MC steps iVgteps and the steps iVthenn disregarded for 
thermalization as used for the determination of the critical Casimir amplitude, for 
lattices with L > 24 and Lii = L/p. Each step corresponds to 1 Metropolis sweep and 
L Wolff single-cluster flips. Runs marked with * have been split into independent runs 
with different random numbers and have then been recollected together: the number 
of steps reported there refers to the cumulated numbers. 



L 


P 


AWs/10 3 


iVthcrm/10 3 


P 


^steps/10 3 


AWrm/10 3 


24 


1/6 


4000* 


200 


1/8 


5100* 


270 


32 


1/6 


5100* 


300 


1/8 


5100* 


270 


40 


1/6 


250 


20 


1/8 


250 


25 


48 


1/6 


1250* 


80 


1/8 


1400* 


140 


G4 


1/6 


2400* 


160 


1/8 


1400* 


80 


24 


1/10 


2000 


100 


1/12 


1750 


20 


32 


1/10 


3280* 


170 


1/12 


1020* 


60 


40 


1/10 


100 


5 


1/12 


90 


20 


48 


1/10 


600* 


30 


1/12 


450 


30 


G4 


1/10 


540* 


44 









Table B2. The number of MC steps N s t ep s, excluded thermalization, used for the 
determination of the order parameter profiles. Each step corresponds to 1 Metropolis 
sweep and L Wolff single-cluster flips. We have sampled the order parameter at 
x/L = ft/8, with k = —16,..., 16, around the chemical step and for z/L — j/8, 
j = 1, ... ,7. The coordinate system is chosen such that the confining surface with 
the spins fixed to form a chemical step corresponds to z — 0, while the homogeneous 
surface with the spins fixed to +1 is located at z — L — 1. In the lateral direction the 
origin x = corresponds, in the lower surface, to spins fixed to +1, the left neighbours 
of which are spins fixed to — 1 . 

L ,9=1/6 p=l/8 ,9=1/10 ,0=1/12 

7V stcps /10 3 iV stcps /10 3 7V stcps /10 3 TVsteps/10 3 



16 


300 


225 


180 


150 


24 


200 


150 


120 


100 


32 


100 


115 


90 


75 


40 


120 


90 


75 


60 


48 


100 


75 


60 


50 



Appendix B. Monte Carlo simulations 

In this appendix we report certain technical details of the Monte Carlo simulations we 
have performed. As explained in section |3j the evaluation of the Casimir force is carried 
out in two steps. First, we determine the thermal average (H2 — "Hi) a which appears in 
eq. ( 1551) . This is done by a standard Monte Carlo simulation for the ensemble given by 
the crossover Hamiltonian Tix defined in eq. (1521) . We implement a combination of the 
standard Metropolis and Wolff cluster algorithms: each MC step consists of 1 Metropolis 
sweep over the entire lattice in lexicographic order and L Wolff single-cluster flips; L 
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denotes the slab thickness including the two surfaces of fixed spins, so that there are L— 2 
layers of fluctuating spins. As random number generator we use the double precision 
SIMD-oriented Fast Mersenne Twister (dSFMT) [97] . Some details of the simulations 
performed at the critical temperature are reported in table IBll Our Simulations have 
been performed by using various clusters. At the critical point they took approximately 
11.5 single-CPU years on a Intel Xeon^ tm ^ E5450 running at 3 Ghz, 11.5 single-CPU 
years on an IBM Power6 ( - tm - ) 575 running at 4.7 Ghz, 3 single-CPU years on a AMD 
Opteron( tm ) 852 running at 2.6 Ghz and 3 single-CPU years on a AMD Opteron^ tm ) 248 
running at 2.2 Ghz. The simulations off the critical point took approximately 38 single- 
CPU years on a Intel Xeon( tm ) E5450 running at 3 Ghz, 47 single-CPU years on an IBM 
Power6( tm ) 575 running at 4.7 Ghz and 6 single-CPU years on a AMD Opteron( tm ) 852 
running at 2.6 Ghz. 

The Metropolis update acts only on the fluctuating spins. Upon implementing the 
Wolff cluster algorithm, particular care has to be taken of the boundary spins, i.e., the 
fluctuating spins located at z — 1 and z = L — 2 which are connected to the two confining 
surfaces, where the spins are either fixed to +1 or to —1. The interactions of the latter 
ones with the surface is given by the surface contribution to the Hamiltonian 

ies + ieS- 

where S + and S*_ are the lattice sites of the boundary spins which are connected to a 
surface spin fixed to +1 and — 1, respectively. The presence of the interaction described 
by eq. (IB.lj) does not allow one to straightforwardly implement the Wolff algorithm. In 
this respect we introduce two fictious spins S p = 1 and S m = — 1 which correspond to 
the spins on the confining surfaces fixed to +1 and to —1, respectively, and we rewrite 
eq. ( IB.ip as 

Us = -j3S p ^2 Si- (3S m Y, Si. (B.2) 

We now want to promote S p and S m to actual fluctuating spins. In order to do so, we 
enlarge our phase space by adding those configurations of spins which correspond to 
reversed signs of S p and S m , while mantaining the relative sign of S p and S m constant. 
In other words, we allow the boundary conditions on the confining surface to be flipped, 
so that the phase space is the union of the phase space corresponding to the geometry 
of figure [1] and the phase space of the "reversed" geometry where the spins at the upper 
surface are fixed to —1 and those at the lower surface are fixed to +1 (x < 0) or to — 1 
(x > 0). The complete partition function of this system reads: 



Z = Yl exp 

{C} 



\{ij) i&s+ ieS- ) i 

fi -> +oo, (B.3) 

where the sum is over the enlarged phase space {C} = {Si = ±l,0}x{S m , S p = ±1} and 
an infinitely strong antiferromagnetic bond between S p and S m ensures that S p = —S m . 
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The partition function given in eq. (1B.3|) is now suitable for applying a cluster routine. 
At the beginning of the simulation we set S p = 1 and S m — — 1. The cluster routine 
selects randomly a spin among the fluctuating ones. If the spin is different from the 
cluster is expanded around such a spin. According to the partition function in eq. (1B.3J) . 



if the cluster reaches one of the confining surfaces, the surface spins are flipped and the 
expansion of the cluster continues for all the spins interacting with the confining surfaces. 
This cluster move only acts on non-zero spins. Ergodicity is obtained by supplementing 
the dynamics with Metropolis sweeps. We mention that a full cluster algorithm is 
possible for the special value D = In 2 [98J . The critical Casimir force resulting from the 
partition function in eq. ( IB. 3 1) is the average of the force for the geometry of figure Q] 
and the "reversed" geometry. Since the force in the two cases is identical, the result is 
indeed the desired critical Casimir force. In fact, upon inspecting eq. ( 149]) one realizes 
that the observable H2 — "Hi, which appears in eq. ( |55j) . is invariant under spin flip. 

The second step of the method consists of the numerical integration in eq. (|55|) . For 
this purpose we employ the Gauss-Kronrod quadrature formula with 15/7 points [99]. 
With this quadrature it is sufficient to sample the integrand in eq. (155]) at 15 points in 
A, the positions of which are fixed. The integral can be estimated using the full set of 
sampled points or using a subset of 7 points. By comparing the results of the numerical 
integration based on 15 and 7 points, respectively, we checked that the systematic 
error due to the discretization of the integral in eq. (155]) is smaller than the statistical 
uncertainty of the Monte Carlo data. 

The order parameter profiles presented in section [6] have been obtained by a 
standard Monte Carlo simulation for the Hamiltonian given in eq. (14*9]) . In table lB2l 
we report some details of these simulations. In this case a cluster flip which involve 
the confining surfaces changes the sign of the magnetization. Accordingly the correct 
profiles are obtained by keeping track of such flips. 
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